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ABSTRACT 

A practical method for identifying linear time in- 
variant systems on the basis of arbitrary input-output 
records is reviewed and extended to handle the case where 
the system is not initially in the zero state. The method 
is implemented using a digital computer program composed 
of a numerical integration subroutine and a subroutine for 
solving overdetermined sets of linear algebraic equations. 
Several examples are presented to demonstrate the accuracy 
and present capabilities of the procedure. 
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I. INTRODUCTION 



A. THE IDENTIFICATION PROBLEM 

In order to apply any of the modern techniques of 
control system design one must first have a mathematical 
model of the system to be controlled. The form of this 
model will depend on the design methods to be employed as 
well as on the physical characteristics of the system. 

Since most of the theory on the analysis and design of 
control systems is based either on the state space or 
transform representation of systems the vast majority of 
mathematical models will consist of either a set of state 
equations or a transfer function. 

Once it has been decided what basic form the mathe~ 
matical model should take the problem of determining the 
numerical values o-f the parameters arises. Parameter 
values can often be determined from the laws of physics 
and the data supplied by a manufacturer or obtained through 
testing. This is not always the case however. Occasionally 
the law's of physics become mathematically intractable or 
are not even applicable. Quite often the values of certain 
key parameters are not available. It is in these cases 
that the identification problem arises. 

A common problem in engineering is that of determining 
the output of a system based on a knowledge of the system 
model, input, and initial conditions. The identification 
problem is similar to this but here the unknown quantity is 
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the system model. The input and output are assumed to be 
known. For the purposes of this thesis the identification 
problem can be stated as follows: 

Given - a record of the input and output of a system 
over some finite period of time, 

Find - a mathematical model and the numerical values 
of the model parameters in such a way that the 
model will accurately describe the behavior 
of the system. 

It should be kept in mind that the problem of identi- 
fying a system solely on the basis of input and output data 
(the so called "black box" problem) is very rare in engi- 
neering. Even in the case where none of the model parame- 
ters are known one will more than likely have a fair idea 
of whether the system is linear or nonlinear, time varying 
or time invariant, what the order of magnitude of the domi- 
nant time constants is, and what types of inputs and out- 
puts are to be expected. For this reason most engineering 
identification problems fall into the "grey box" category. 
This distinction may seem trivial but virtually all iden- 
tification techniques rely heavily on knowledge of the 
characteristics and quantities mentioned above. 

Although the control systems literature on system 
identification is quite vast there are no known identifica- 
tion schemes capable of handling all identification prob- 
lems effectively. Choosing a method suitable for a given 
problem can become a formidable task. A paper summarizing 
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most of the common approaches to the problem of identifying 
lumped parameter systems has been published by Nieman, 

Fisher, and Seborg [1], A good discussion of the indus- 
trial applications of various methods has been published 
by Eykhoff, et al . [2]. 

One common approach used in linear system identifica- 
tion is that of obtaining the impulse response of the sys- 
tem. Mishkin and Haddad [3] have developed a technique 
for finding the impulse response based on samples of the 
system output and its derivatives. A technique for esti- 
mating the impulse response on the basis of noisy input and 
output samples has been developed by Levin [4] and Kerr and 
Surber [5]. Turin [6] and Lichtenberger [7,8] have used a 
matched filter to obtain an identification. The use of a 
white noise or binary test signal and crosscorrelation has 
been suggested by Hill and McMurtry [9] . The noise limita- 
tions of the sample approximation, matched filter, and 
crosscorrelation identification techniques have been in- 
vestigated by Lindenlaub and Cooper [10]. 

Another common approach to the identification problem 
is to determine the coefficients of the differential or 
difference equation which describes the system. Kumar and 
Sridhar [11] have employed the method of quasilinearization 
with some success. -Nagumo and Noda [12] have developed a 
learning approach to the problem. ' Bass [13] has developed 
a method which uses modulating functions and works well in 
the presence of noise. Astrom and Bohlin [14] have developed 
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a statistically optimal method of determining differential 
equation parameters known as the "maximum likelihood method." 
A similar method which is not optimal but is considerably 
simpler computationally has been developed by Peterka and 
Smuk [15,16]. It is known as the "prior knowledge fitting 
method." An algorithm for determining state variable models 
of sampled data systems has been proposed by Ho and Kalman 
[17]. The algorithm performs quite well in the presence 
of noise, and has been extended to continuously operating 
systems by Eldem [18] . 

Methods of identifying nonlinear and distributed param- 
eter systems are usually limited to specific types of sys- 
tems or to specific types of nonlinearities. This is 
undoubtedly due to the wide variety of nonlincarities en- 
countered in physical systems and the difficulty of finding 
a model capable of characterizing them all. Shinbrot [19], 
Mowery [20], Fairman [21], and Bellman, et al. [22] have 
all approached the problem of identifying nonlinear systems 
by assuming a particular form of differential equation is 
capable of describing the system and then developing methods 
around the form of differential equation chosen. Another 
common approach to nonlinear system identification is that 
of representing a system by a suitable functional polyno- 
mial relating the input and output. Hsieh [23] uses this 
approach and a steepest descent algorithm to solve the iden- 
tification problem. Similar approaches have been taken by 
Simpson [24], Bose [25,26], and Hubbell [27]. 
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Identification methods vary widely with respect to how 
much must be known about the system before the method can 
be applied. Some identification techniques require that 
prior estimates of all system parameters be available. 

Many methods restrict the allowable system inputs to a 
family of testing functions such as steps or binary pulses. 
In general, the less that is known about a given system and 
the tighter the constraints on the kind of signals which 
may be applied as inputs the more difficult it is to find 
a method capable of accomplishing the identification. 

B. OBJECTIVES OF THIS PAPER 

This paper will present a study of an identification 
technique originally suggested by Diamessis [28]. It is 
designed to identify lumped linear time invariant systems 
but has been extended by Diamessis [29] and Wang [30] to 
handle certain types of nonlinearities. The technique re- 
quires a knowledge of the system input and output over some 
finite time interval as well as a rough estimate of the sys- 
tem order. The system input need not be restricted to a 
class of testing functions, it can be completely arbitrary. 

Unlike some identification techniques v^rhich require 
the calculation of derivatives of the input and output, the 
technique to be presented requires only integrals of the 
input and output. The advantages of numerical integration 
over differentiation are well known. Since any zero mean 
noise component on the input or the output tends to be 
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attenuated greatly by the integration process the system 
identification can be more accurate than the raw data used 
to accomplish it. • 

The remainder of this thesis is divided into three 
major sections. In Chapter II the theoretical development 
of the identification technique is given. A method for 
identifying the initial conditions of the unknown system 
is also presented. Chapter III presents a method for im- 
plementing the techniques developed in Chapter II. Par- 
ticular attention is given to the choice of numerical 
methods and to efficient programming techniques. Several 
examples are presented to demonstrate the capabilities of 
the method. In Chapter IV several recommendations for 
further study are made. Conclusions concerning the accu- 
racy and present limitations of the technique under con- 
sideration are also discussed. Following the conclusions 
a complete listing of all computer programs used in the 
thesis is given. 
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II. IDENTIFICATION BY MULTIPLE INTEGRATIONS 



A. GENERAL APPROACH 

The development which follows is similar to the devel- 
opment given by Diamessis [28] in 1965. There are a few 
notable differences however. The development given by 
Diamessis is restricted to the case where all initial con- 
ditions are zero. This is a rather serious restriction 
since it may be difficult or impossible to find a point 
where the system is in the zero state if the systems opera- 
tion is not to be disturbed. Zero initial conditions will 
not be assumed in the development which follows. A method 
for solving for the unknown initial conditions will be 
presented. Diamessis proposed the formulation of a uniquely 
determined set of linear algebraic equations with the model 
parameters as unknowns. This development will make use of 
overdetermined sets of linear algebraic equations with the 
model parameters and initial conditions as unknowns. The 
overdetermined set of equations will then be solved using 
the method of least squares. It will be sho\m that this 
results in a more accurate identification when the accuracy 
of the available data is limited and the order of the sys- 
tem is unknown. 

Any single input, single output, lumped parameter, 
linear, time-invariant system can be described by a linear 
ordinary differential equation with constant coefficients. 
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The basic form of this equation is given in Equation (1) 
along with a set of initial conditions. 
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where ; 

u(t) = system input 
y(t) = system output 

The identification problem to be treated here consists 
of determining n, m, and the various coefficients of the 
differential equation on the basis of input and output 
records taken over some arbitrary time interval. The ini- 
tial conditions will be assumed to be unknown. 

Taking the Laplace transform of Equation (1) yields 
Equation (2) . 
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The coefficients account for the contributions of the 
initial conditions. Dividing Equation (2) by s^^^ is 
equivalent to integrating n+1 times in the time domain. 
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Taking the inverse transform of Equation (3) results in 
Equation (4) . 
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Since the system is time invariant nothing has been lost by 
setting the lower limit on the integrals equal to zero. 
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Rearranging terms in Equation (4) and placing all terms which 
depend on u(t) , y(t) , or t in brackets results in Equation 
(5). 

(5) 
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Since records o£ the input and output are assumed to 
be known a linear algebraic equation with the system param- 
eters and initial condition terms as unknowns can be formu- 
lated by performing the indicated multiple integrations 
from zero to some time tj,. A set of 2n+m+l equations can 
be obtained by letting tj^ take on 2n+m+l different values. 
Assuming that the equations are linearly independent it 
will now be possible to solve for the n+m+1 differential 
equation coefficients and the n initial condition terms. 

So far nothing has been said of how n and m are de- 
termined. Theoretically it should be possible to use any 
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n' and m' greater than the actual order of the system under 
study. If the order of the system is n with m input coef- 
ficients then one would expect the following: 

a^ = 0 for ( 0 < i < n'-n-l ) 

= 0 for ( 0 < i < m’-m-l ) 

with n'>n and m'>m. 

The model should essentially reduce itself to the right order 
by setting nonessential terms equal to zero. 

Unfortunately the situation is not quite this simple. 

Due to the finite accuracy of all experimental data the lin- 
ear equations will not have an exact solution. For certain 
types of inputs the linear equations will not even be lin- 
early independent. These problems can be overcome to some 
extent by formulating more than 2n+m+l equations which are 
required. The overdetermined set can then be solved using 
the method of least squares. If the limited accuracy of 
the experimental data can be attributed to round off errors 
then integrating the data over a time interval much greater 
than the sampling interval should remove much of the uncer- 
tainty in the linear equation coefficients. This is be- 
cause the roundoff process can be modeled as a zero mean 
white noise process. The integral of the noise will ap- 
proach zero as the period of integration increases. 

Even the measures mentioned above will not solve the 
problem completely however. Due to the finite precision 
used to represent numbers in a computer and the iterative 
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nature of the numerical methods used to solve overdeter- 



mined sets of linear equations it is impossible to obtain 
zero as a solution for any parameter. If the parameter 
should be zero the numerical method will return a very 
small but nonzero value. Although this will result in an 
incorrect estimate of the system order the error will not 
be serious in most cases. This is because the small co- 
efficients of the terms which should nonexistent will make 
their effect negligible. Examples in Chapter III will 
demonstrate this point. 

Although the development in this section has assumed 
a differential equation model of the system the same results 
could have been obtained if a transfer function or state 
variable model had been assumed. A simple, rearrangement of 
terms in Equation (2) results in Equation (6) . 



Y(s) b„s"^ + + b,s + b^ (6) 

m 1 o ^ 



U(s) s^ + a ,s^ ^ + + a, s + a 
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By defining a few new terms. Equation (7) results. 

Y(s) K ( s''^ + c , s’^ ^ + + c. s + c (7) 

^ ^ m-1 1 o ^ ^ 



U(s) s^ + a ,s*^ ^ + + a,s + a 

n-1 1 o 

where K = b 

m 

c. = b. / b for 0 < i < m. 

1 1 m - - 



A set of state equations may be formulated in a similar 
fashion. One convenient state variable representation is 
given in Equation (8) . It is based on the phase variable 
form of system representation. 
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Note the simple correspondence between the terms in the trans- 
fer function and the terms in the state equations. 



B. IDENTIFYING INITIAL CONDITIONS 

In many identification problems it is desirable to com- 
pare...the system model with the actual system by exciting the 
system model with the same input data that was used in the 
identification. By comparing the output of the model with 
the output of the system a rough idea of the accuracy of the 
model can be obtained. This will not be possible however 
unless the initial conditions at the beginning of the input- 
output record are all known. 

From Equation (5) it can be seen that when the linear 
algebraic equations are solved for the unknown model parameters 



17 



the initial condition terms are also found. By taking 
the Laplace Transform of Equation (1) and specifically 
writing in the contributions of the individual initial 
conditions a relationship between the g^ terms and the ini- 
tial conditions can be found. 
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Unfortunately Equation (9) requires the knowledge of 
m-1 derivatives of the input function u. It may be neces- 
sary to calculate m-1 derivatives of the input using numer- 
ical techniques. This could cause the model and system 
output to differ slightly at the beginning of the output 
record but as the natural response dies out the records 
should converge. It may be possible to avoid this diffi- 
culty in many cases by choosing the beginning of the input- 
output record at a point where the input is relatively 
constant . 

C. SIMPLIFICATIONS WITH ZERO INITIAL CONDITIONS 

In many problems it will be possible to exercise com- 
plete control over the input to the system under study. If 
it is possible to obtain an input-output record beginning 
when the system is in the zero state it will be possible to 
simplify the identification procedure. Since the g^ terms 
will all be zero if the system is in the zero state they 
' need not be included in the formulation of the linear alge- 
braic equations. This will reduce the number of unknowns 
from 2n+m+l to n+m+1. 

It should be noted that additional information can 
often be incorporated into the identification procedure in 
order to simplify the problem. For example if the steady 
state gain constant were known the number of unknowns could 
be reduced by one. It is usually a simple matter to tell 
whether a system has poles or zeroes at the origin. This 
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piece o£ information can easily be used to simplify the 
identification procedure still further. As a general rule, 
the fewer the unknowns the more accurate the identification. 
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III. IMPLEMENTATION 



A. NUMERICAL METHODS 

The identification technique presented in Chapter II 
can be broken into two basic steps. The first step con- 
sists of performing the multiple integrations of the input 
and output and forming the overdetermined set of linear 
algebraic equations. The second step consists of solving 
the linear equations for the unknown model parameters. It 
is a distinct advantage of the identification technique 
under study that each of these steps can be carried out by 
subprograms that are readily available in virtually all 
modern computer centers. 

Step one can be handled by any numerical integration 
subroutine. Although there are quite a few highly sophis- 
ticated numerical integration procedures available, trape- 
zoidal integration will give better results in most 
applications. There are several reasons why this is true. 
First of all, most of the more complex integration tech- 
niques perform poorly when the function being integrated 
is discontinuous. Since control system inputs are often 
discontinuous and since such discontinuities are quite 
desirable from an identification standpoint, complex inte- 
gration techniques are usually undesirable. Even \>?hen the 
functions to be integrated are continuous the slight im- 
provement in accuracy offered by more advanced methods is 
not enough to justify the tremendous increase in computation 
al load associated with their use. 
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Step two, the solution of the set of overde termined 
linear equations, is a classical problem in several fields 
of mathematics and engineering. Unfortunately most of the 
classical techniques for solving such problems are not 
practical. They tend to magnify the errors introduced by 
the finite precision of the computer to the point where 
the solution is meaningless. Fortunately several modern 
methods are available that display more acceptable behavior. 
The method used in this paper was developed by Golub [31] 
in 1965. The basic approach is to triangularize the coef- 
ficient matrix by performing a Choleski decomposition. The 
decomposition is accomplished by applying repeated House- 
holder transformations [32] . Once the coefficient matrix 
has been triangulari zed the unknowns can be obtained by 
back substitution. The method is quite stable numerically 
and is capable of handling illconditioned coefficient ma- 
trices . 

B. CHARACTERISTICS OF GOOD INPUT-OUTPUT RECORDS 

The accuracy with which a system can be identified is 
strongly dependent on the input-output record used in the 
identification. Since parameters are identified on the 
basis of their effect on the output it will be impossible 
to identify a parameter unless its effect is measureable. 

If the input to a system has a frequency spectrum that is 
more or less uniform over the frequency range of interest 
the identification will probably be very good. If the 



22 



frequency spectrum of the input is confined to a narrow 
band of frequencies the identification will probably be 
very bad. It is well known that signals with sharp dis- 
continuities have a broader bandwidth than slowly varying 
continuous signals. For this reason input-output records dis- 
playing discontinuities and rapid time variations should 
be chosen. 

If step or ramp inputs are used in the identification 
the value of the initial conditions will have to be known 
and incorporated into the set of linear equations. Since 
the initial conditions will usually be zero when these 
types of inputs are used this should not cause any diffi- 
culties. The reason for this difficulty lies in the nature 
of the initial condition coefficient terms. Tlie integral 
coefficients of these terms are steps, ramps, and higher 
order polynomials in time. If the input is a step or a 
ramp the coefficients of several model parameters will also 
be steps, ramps, and higher order polynomials in time. 

There will therefore be a direct relationship between model 
parameter coefficients and initial condition term coeffi- 
cients. This will result in the linear equations having an 
infinite number of solutions due to the linear dependence 
between all the individual equations in the set. Step and 
ramp inputs must therefore be avoided when the system ini- 
tial conditions are unknown. 

Since analog system data will have to be converted to 
digital form a suitable sampling interval will have to be 
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chosen. Experimental results have shown that a sampling 
rate ten to one hundred times shorter than the shortest 
system time constant works quite well. Lower sampling 
rates may introduce inaccuracies in the location o£ high 
frequency poles and zeroes. 

C. PROGRAM DESCRIPTION 

In order to evaluate experimentally the characteristics 
of the identification procedure under study a set of digital 
computer programs was developed. The main identification 
program is a direct implementation of the procedure develop- 
ed in Chapter II. Subroutine SYSTEM is a simulation program 
that was written to generate input-output data for the main 
program to process. 

In the beginning of the identification program several 
important parameters are defined. NP and NZ are rough es- 
timates of the number of poles and the number of zeroes in 
the system to be identified. KPTMAX is the number of sam- 
ple points available in the input-output record. Each sam- 
ple point will consist of the time (T) , the input amplitude 
(R) , and the output amplitude (C) . IPTS is the number of 
sample points that will separate successive linear equa- 
tions. In other words, every time the total number of points 
read in is a multiple of IPTS a linear equation will be gen- 
erated. The total number of linear equations that will be 
generated is equal to KPTMAX/IPTS. When all of the linear 
equations have been formed subroutine DLLSQ is called. This 
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subroutine is an implementation of the Golub algorithm for 
solving overdetermined sets of linear equations. 

Subroutine DLLSQ returns the values of the system 
model parameters and initial condition parameters of Equa- 
tion (2) . In order to find the poles and zeroes of the 
system the output of DLLSQ is fed into RTPLSB. RTPLSB is 
a polynomial root finder which .uses a combination of the 
Newton-Raphson and Bairstow methods to find the poles 
and zeroes of the system. 

The remainder of the main identification program is 
devoted to output. The results of the identification are 
given in both transfer function and state variable form. 

The state variable form is referenced to the format used 
in Equation (8) . 

Subroutine SYSTEM reads in the transfer function of a 
system and computes the system output based on a set of 
arbitrary initial conditions and an arbitrary input wave- 
form. Each time subroutine SYSTEM is called by the iden- 
tification program it returns three numbers, the time T, 
the input waveform amplitude R, and the output waveform 
amplitude C. Subroutine SYSTEM prints out the transfer 
function and state variable representation of the system 
it is simulating so that the accuracy of the identification 
can be determined. 

Input-output recrods obtained from physical systems 
are rarely accurate to more than three or four significant 
digits. The data generated by subroutine SYSTEM is therefore 
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rounded off by subroutine ROUND before being passed to the 
identification program. The number of significant digits 
in the data returned to the identification program may be 
varied by changing the value of the parameter NA in the 
simulation subroutine. 

D. EXAMPLES 

The following examples demonstrate ho;^^ the accuracy of 
an identification depends on the accuracy of the input- 
output record, the sampling period, and the input waveform. 
They also show how the order of the system can be deter- 
mined from a trial identification run using an estimated 
order greater than the actual order of the system. 

Example one illustrates the relationship between the 
accuracy of the input-output record data and the resulting 
identification. In order to minimize the effect of other 
factors all initial conditions were set equal to zero, a 
step input was used, and n' and m' were set equal to n and 
m. Example one demonstrates the fact that there is a di- 
rect (almost linear) relationship between the accuracy of 
input-output data and the accuracy of the identification. 
Note that even when the input-output record contained only 
two significant digits the identification of the system 
poles and zeroes was within about three percent of their 
exact values . 

Example two illustrates the relationship between the 
sampling period used with the input-output records and the 
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the accuracy o£ the identification, input records contain- 
ing discontinuities or rapid time variations should be 
chosen . 

Examples four through ten demonstrate what happens 
when n’ and m’ are greater than n and m. In each example 
an identification is performed using an n' and m' greater 
than n and m. Using the information obtained from this 
identification a new value of n' and m' is determined. 

These new estimates are then used to perform a second iden- 
tification . 

In examples four and five the error in the estimate of 
m' was made larger than the error in the estimate of n'. As 
a result of the relative values of these two errors all the 
excess poles cancelled with excess zeroes but since there 
were more excess zeroes than excess poles some excess zeroes 
remained. Note however that for frequencies lower than the 
sampling rate the excess zeroes have little or no effect 
on the behavior of the identified system. Experiments have 
shown that excess zeroes that do not cancel with excess 
poles are always of a frequency comparable to or higher 
than the sampling rate. Using this principle and estima- 
ting new values for n' and m' results in an identification 
which has the correct number of poles and zeroes and is 
more accurate than the first identification. 

In examples six, seven, and eight the error in the 
estimate of n' was equal to the error in the estimate of 
m'. As a result, all excess poles and zeroes cancelled 
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with each other leaving a system o£ the correct order. 

Note that by repeating the identification with the correct 
value of n' and m’ it was possible to improve the accuracy 
of the identification. 

In examples nine and ten the error in the estimate 
of n' was made greater than the error in the estimate of 
m' . As a result, all excess zeroes cancelled with excess 
poles but some excess poles remained. Note that for fre- 
quencies lower than the sampling rate the excess poles 
have negligible effect. Experiments have shown that excess 
poles that do not cancel with excess zeroes are almost al- 
ways of a frequency comparable to or greater than the sam- 
pling rate. Using this principle and estimating new values 
for n' and m' resulted in a correct identification of the 
simulated systems. 

Experiments have shown that identifications involving 
uncancelled excess zeroes are generally more accurate than 
identifications involving uncancelled excess poles. For 
this reason it is best to set m’ close to n' when identi- 
fying an unknown system. 

Using the experimental findings described above a 
simple procedure for determining n and m can be formulated. 
First, guess an n’ which is greater than the order of the 
system to be identified. This should not be too difficult 
in most engineering identification problems. Let m' be 
equal to n’ or m’-l. This will guarantee complete cancel- 
lation of all excess poles and zeroes or partial cancellation 
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of excess poles and zeroes with excess zeroes remaining. 
After making a trial identification with the values of n' 
and m' mentioned above, estimate new values for n' and m' 
based on the reasoning in the examples. The neiv values 
of n' and m' should now be correct. By performing the 
identification with these values of n’ and m' a good iden- 
tification should result. 
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EXAMPLE 1 



SYSTEM TO BE IDENTIFIED 



SYSTEM TRANSFER FUNCTION 

POLES REAL 

1 - 1.0000000000 

2 - 10.0000000000 

3 -lOO.COOOOOOOOO 



IMAGINARY 

0.0 

0.0 

0.0 



J 

J 

J 



ZEROES REAL 

1 -3.0000000000 

2 -30.0000C000C0 



IMAGINARY 
0.0 J 

0.0 J 



GAIN CONSTANT = 10.0000000000 



SYSTEM state VARIABLES (PHASE FORM) 



A 


VECTOR 






1 


ICOO.OOOOOCOOOO 




2 


lllC.COOOOOOOCO 




3 


lll.OCOOOOOOOO 


B 


VECTOR 






3 


10.0000000000 


C 


VECTOR 






1 


90.0000000000 




2 


33.C0000000C0 




3 


1 . cooocoooco 
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I 



I 





IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES REAL 

1 -100.0001690225 

2 -0.9999970711 

3 -9.999988A241 



ZEROES REAL 

1 -2.9999937654 

2 -30.0000138488 



IMAGINARY 
0.0 J 

0.0 J 

0.0 J 

IMAGINARY 
0.0 J 

0.0 J 



GAIN CONSTANT = 10.0000038147 



SYSTEM STATE VARIABLES (PHASE FORM) 



A 


VECTOR 






1 


999.9976037596 




2 


1110.0G03679024 




3 


111.0001545177 


B 


VECTOR 






3 


10.0000038147 


C 


VECTOR 






1 


89.9998545091 




2 


33.0000076143 




3 


1.0000000000 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



0.60000 SEC 
0.20000 MSC 
120 

8 DECIMAL PLACES 
UNIT STEP 
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r 



•W' 




IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 




1 


-100.0000291125 


o 

• 

o 


J 


2 


-1.0000011522 


o 

• 

o 


J 


3 


-10.0000016550 


o 

• 

o 


J 


ZEROES 


REAL 


IMAGINARY 




1 


-3.0000028353 


o 

• 

o 


J 


2 


-30.0000250480 


o 

• 

o 


J 



GAIN CONSTANT = 9.99999332A3 



SYSTEM STATE VARIABLES (PHASE FORM) 
A VECTOR 

1 1C00.0016C88045 

2 1110.0006141285 

3 111.0000319197 

B VECTOR 

3 9.9999933243 

C VECTOR 

1 90.00C1602028 

2 33.0000278833 

3 l.COOOOCOCOO 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



0.60000 SEC 
C. 20000 MSC 
120 

7 DECIMAL PLACES 
UNIT STEP 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES REAL 

1 -1C0.000C683416 

2 -1.0000118079 

3 -10.0000329083 



ZEROES REAL 

1 -3.0000272949 

2 -30.0000668753 



IMAGINARY 
0.0 J 

0.0 J 

0.0 J 

IMAGINARY 
0.0 J 

0.0 J 



GAIN constant = 9.9999847412 



SYSTEM STATE VARIABLES (PHASE FORM) 

A VECTOR 
1 
2 
3 

B VECTOR 

3 9.9999847412 

C VECTOR 

1 90.0010194754 

2 33.0000941702 

3 l.OOOOOCOOOO 



lOCO. 0157822215 
1110.0053743721 
111.0001130579 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



0.60000 SEC 
0.20000 MSC 
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120 

6 DECIMAL PLACES 
UNIT STEP 



IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 




1 


-99.9968491528 


0.0 


J 


2 


-1.C000017165 


o 

• 

o 


J 


3 


-9.9998390264 


o 

• 

o 


J 


ZEROES 


REAL 


IMAGINARY 




1 


-2.9999754422 


o 

• 

o 


J 


2 


-29.9992863605 


o 

• 

o 


J 



GAIN CONSTANT = 9.9998254776 



SYSTEM STATE VARIABLES (PHASE FORM) 

A VECTOR 
1 
2 
3 

B VECTOR 

3 9.9998254776 

C VECTOR 

1 89.9971223661 

2 32.9992618027 

3 1.0000000000 



999.9541111369 

1109.9492716708 

110.9966898958 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



0.60000 SEC 
0.20000 MSC 
120 

5 DECIMAL PLACES 
UNIT STEP 
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I 






IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 




1 


-99,9927835443 


o 

• 

o 


J 


2 


-0.9991304951 


o 

• 

o 


J 


3 


-9,9980805636 


o 

• 

o 


J 


ZEROES 


REAL 


IMAGINARY 




1 


-2. 9981092090 


o 

• 

o 


J 


2 


-29. 9981888298 


o 

• 

o 


J 



GAIN CONSTANT = 9,9994869232 



SYSTEM STATE VARIABLES (PHASE FORM) 

A VECTOR 
1 
2 

3 - 

B VECTOR 

3 9,9994869232 

C VECTOR 

1 89,9378461845 

2 32,9962980388 

3 1.0000000000 



998. 8666303721 
1109.6311321656 
110,9899946030 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT function 



C.6U000 SEC 
0.20000 MSC 
120 

4 DECIMAL PLACES 
UNIT STEP 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 




1 


-99.9937910594 


o 

• 

o 


J 


2 


-0.9996622783 


o 

• 

o 


J 


3 


-10.0015833055 


0.0 


J 


ZEROES 


REAL 


IMAGINARY 




1 


-2.9996507024 


o 

• 

o 


J 


2 


-30.0061026888 


o 

• 

o 


J 



GAIN CONSTANT = 9.9986057281 



SYSTEM STATE VARIABLES (PHASE FORM) 

A VECTOR 
1 
2 
3 

B VECTOR 

3 9.9986057281 

C VECTOR 

1 90.0078270058 

2 33.CC57533912 

3 l.OCOOOGOOGO 



999.7584771470 

1110.0544578539 

110.9950366432 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



0.60000 SEC 
0.20000 MSC 
120 

3 DECIMAL PLACES 
UNIT STEP 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL IMAGINARY 


1 


-98. 8C04578528 0.0 J 


2 


-0.9763637189 0.0 J 


3 


-9.8862584052 0.0 J 


ZEROES 


REAL IMAGINARY 


1 


-2.9406610975 0.0 J 


2 


-29.6497245014 0.0 J 



GAIN CONSTANT = 9.945189476C 



SYSTEM STATE 


VARIABLES (PHASE FORM) 


A VECTOR 
1 


953.6797208413 


2 


1082. 8846233636 


3 


109. 6630799769 


B VECTOR 
3 


9.9451894760 


C VECTOR 
1 


87.1897913926 


2 


32. 5903855989 


3 


1.0000000000 



PROGRAM PARAMETERS 

RECORD LEGNTH = 0.60000 SEC 

SAMPLING PERIOD = 0.20000 MSC 

EQUATIONS FORMED = 120 

PRECISION OF DATA = 2 DECIMAL PLACES 

INPUT FUNCTION = UNIT STEP 
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EXAMPLE 2 



SYSTEM TO BE IDENTIFIED 



SYSTEM TRANSFER FUNCTION 

POLES REAL 

1 - 1.0000000000 

2 - 10.0000000000 

3 -100.0000000000 



IMAGINARY 



0.0 J 

0.0 J 

0.0 J 



ZEROES REAL 

1 -3.O000C00OC0 

2 -30.0000000000 



IMAGINARY 
0.0 J 

0.0 J 



GAIN CONSTANT = 10. OOODCOOOOO 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 lOOO.OOOOOOOOCO 

2 lllO.COOOOOOOOC 

3 111.0000000000 

B VECTOR 

3 10.0000000000 



C VECTOR 

1 90.COCOCOOOOO 

2 33.0000C00000 

3 1.0000000000 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 

1 

2 

3 



REAL 

-101.0610717819 

-1.0000331312 

-10.0017011716 



ZEROES REAL 

1 -3.0002183624 

2 -30.0553440030 



imaginary 

0.0 J 

0.0 J 

0.0 J 

IMAGINARY 
0.0 J 

0.0 J 



GAIN CONSTANT = 10.0887174606 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 1010,8161234948 

2 1121.8490926385 

3 112.0628060847 

B VECTOR 

3 10.0887174606 

C VECTOR 

1 90.1725949663 

2 33.0555623654 

3 1.0000000000 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OP DATA 
INPUT FUNCTION 



0.60000 SEC 
0.10000 MSC 
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120 

3 DECIMAL PLACES 
UNIT STEP 



IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 

1 -105 

2 -1 

3 -10 

ZEROES 

1 -3 

2 -30 

GAIN CONSTANT = 



REAL 

3244028971 

0000069733 

0258525504 

REAL 

0005906385 

4710561855 

10.3944911957 



IMAGINARY 
0.0 J 

0.0 J 

0.0 J 

IMAGINARY 
0.0 J 

0.0 J 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 1055.9742969792 

2 1171.3179932211 

3 116.3502624208 

B VECTOR 

3 10.3944911957 

C VECTOR 

1 91.4311659362 

2 33.4716468241 

3 l.OOOOOOCOCO 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



1.20000 SEC 
0.200C0 MSC 
120 

3 DECIMAL PLACES 
UNIT STEP 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 




1 


-107. 6394221979 


o 

• 

o 


J 


2 


-1.0000162038 


o 

• 

o 


J 


3 


-10.0077785631 


o 

• 

o 


J 


ZEROES 


REAL 


IMAGINARY 




1 


-3.0001692900 


o 

• 

o 


J 


2 


-30.3047449110 


o 

• 

o 


J 



GAIN CONSTANT = 10.6635513306 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 1C77. 2489572831 

2 1194.8806091150 

3 118.6472169649 

B VECTOR 

3 10.66355133C6 



C VECTOR 

1 90.9193650217 

2 33.3C49142CC9 

3 l.COCOOOOOOO 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



= 3.00000 SEC 
= 0.50000 MSC 
= 120 

= 3 DECIMAL PLACES 
= UNIT STEP 
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EXAMPLE 3 



SYSTEM TO BE IDENTIFIED 
SYSTEM TRANSFER FUNCTION 

J 
J 
J 

J 
J 

GAIN CONSTANT = 10.0000000000 



POLES 


REAL 


IMAGINARY 


1 


-1.0000000000 


o 

• 

o 


2 


-10.0000000000 


o 

• 

o 


3 


-100.0000000000 


o 

• 

o 


ZEROES 


REAL 


IMAGINARY 


1 


-3.0000000000 


o 

• 

o 


2 


-30.0000000000 


o 

• 

o 



SYSTFM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 lOOO.COOCOOOOCO 

2 lllO.OOOOOOOUOC 

3 111.0000000000 



B VECTOR 

3 10.0000000000 



C VECTOR 

1 90.0000000000 

2 33.0COOOOOOOO 

3 l.OOOOOOOOCO 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES REAL 

1 -99.6074714180 

2 -1.0627466213 

3 -10.0526991957 



ZEROES REAL 

1 -3.1109954642 

2 -29.8941214051 



IMAGINARY 
0.0 J 

0.0 J 

0.0 J 

IMAGINARY 
0.0 J 

0.0 J 



GAIN CONSTANT = 9.9883327484 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 1064.1536423652 

2 1117.8649236164 

3 110.7229172350 

B VECTOR 

3 9.9883327484 



C VECTOR 

1 93.0004760964 

2 33.0051168692 

3 1.0000000000 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA = 
INPUT FUNCTION 



0.60000 SEC 
C. 10000 MSC 
120 

3 DECIMAL PLACES 
GAUSSIAN- SIGMA=0.010 SEC 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 




1 


-98.9033657608 


o 

• 

o 


J 


2 


-0.9389787239 


o 

• 

o 


J 


3 


-9.8317762550 


o 

• 

o 


J 


ZEROES 


REAL 


IMAGINARY 




1 


-2. 8677895731 


o 

• 

o 


J 


2 


-29.6513565567 


o 

• 

o 


J 



GAIN CONSTANT = 9.9508285522 



SYSTEM STATE VARIABLES (PHASE FORM) 

A VECTOR 
1 
2 
3 

B VECTOR 

3 9,9508285522 

C VECTOR 

1 85.0338511612 

2 32.5191461297 

3 1.0000000000 



913.0589327259 

1074.4957479197 

109.6741207397 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



0.60000 SEC 
C.IOOGO MSC 
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120 

3 DECIMAL PLACES 
GAUSSIAN- SIGMA=0.030 SEC 



IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 




1 


-99, 1859265460 


o 

• 

o 


J 


2 


-0.7535004508 


o 

• 

o 


J 


3 


-9. 5204582353 


o 

• 

o 


J 


ZEROES 


REAL 


IMAGINARY 




1 


-2.5186589038 


o 

• 

o 


J 


2 


-29.4122032416 


o 

• 

o 


J 



GAIN CONSTANT = 9.9813451767 



SYSTEM STATE VARIABLES (PHASE FORM) 

A VECTOR 
1 
2 
3 

B VECTOR 

3 9.9813451767 

C VECTOR 

1 74.0793075760 

2 31.9308621454 

3 l.OOOOOCCOOO 



711. 5270631979 
1026.2057811417 
109.4598852321 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



0.6GC00 SEC 
O.IOOCO MSC 
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120 

3 DECIMAL PLACES 
GAUSSIAN- SIGMA=0.050 SEC 



EXAMPLE 4 



SYSTEM TO BE IDENTIFIED 
SYSTEM TRANSFER FUNCTION 



POLES 


REAL IMAGINARY 


1 


-5.0000000000 0.0 J 


2 


-20.00C00C00G0 0.0 J 


ZEROES 


REAL IMAGINARY 


GAIN CONSTANT 


= 300.0CCOOOOOOO 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 




1 


loo.oocooocoon 


2 


25.0000000000 


B VECTOR 




2 


300 . ooocoocooo 


C VECTOR 




1 


l.OOOOOCCOOC 



INITIAL STATE VECTOR 

1 l.COOOOOOOOO 

2 0.0 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 


1 


-20. 1020603504 


0.0 J 


2 


-0.4280895557 


-6.7540893481 J 


3 


-0.4280895557 


6.7540893481 J 


4 


1.0489471395 


0.0 J 


5 


-4.9974142741 


0.0 J 


ZEROES 


REAL 


IMAGINARY 


1 


211 68. 5357639489 


0.0 J 


2 


-0.4631078779 


-6.7682935708 J 


3 


-0.4681078779 


6.7682935708 J 


4 


1.0470216652 


0.0 J 


GAIN CONSTANT 


= -0.0141926892 





SYSTEM STATE 


VARIABLES <PHASE 


FORM) 


A VECTOR 


1 


-4826. 300213 




2 


3305.020433 




3 


1059.631159 




4 


140.5228441 




5 


24.90670660 




B VECTOR 


5 


-0. 1419268921E 


-01 


C VECTOR 


1 


102C181.159 




2 


-953662.8857 




3 


2390. 647541 




4 


-21168. 64657 




5 


1 . OCOOOOOOO 
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I 





IDENTIFICATION OF UNKNOWN SYSTEM 



INITIAL 

1 

2 

3 

4 

5 



CONDITION (G) VECTOR 
0. 9963029652 
25.01148179 
37. 55868385 
1 11 8.728299 
-1 523. 539669 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



3.00000 SEC 
0.50000 MSC 
120 

3 DECIMAL PLACES 
PIECEWISE CONSTANT 
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IDENTIFICATION OF UNKNOWN SYSTEM 

SYSTEM TRANSFER FUNCTION 

POLES REAL 

1 -4.9994757762 0. 

2 -20.0118644519 0. 

ZEROES REAL 

GAIN CONSTANT = 300.1457519531 



SYSTEM STATE VARIABLES {PHASE FORM) 



A VECTOR 
1 
2 

B VECTOR 
2 

C VECTOR 
1 



100.0488316 

25.01134023 



300.1457520 



1.000000000 



INITIAL CONDITION (G) VECTOR 

1 1.001249200 

2 25.C0380377 



IMAGINARY 
0 J 

0 J 

IMAGINARY 
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EXAMPLE 



5 



SYSTEM TO BE IDENTIFIED 



SYSTEM TRANSFER FUNCTION 

POLES REAL 

1 -2.O00OO0C000 

2 - 2.0000000000 

3 -20.0000000000 



IMAGINARY 
3.0000000000 J 
-3.0000000000 J 
0.0 J 



ZEROES 

1 



REAL 

- 8.0000000000 



IMAGINARY 
0.0 J 



GAIN CONSTANT = 160.0000000000 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 
1 
2 
3 



260.0000000000 
93. OOOOOOOOCO 
24.0000000000 



B VECTOR 
3 



160.C000000000 



C VECTOR 
1 
2 



8.000000C000 

l.OOOOOOOOOOi 



INITIAL STATE VECTOR 

1 2.0000000000 

2 1.0000000000 

3 0.0 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER 

POLES 

1 

2 

3 

4 

5 



FUNCTION 

REAL 

-20.3839954707 

1.0012938647 

1.0012938647 

-1.9998854338 

-1.9998854338 



ZEROES 

1 

2 

3 

4 



REAL 

-7. 6038902797 
447. 1906330529 
1.2922226696 
1. 2922226696 



IMAGINARY 
0.0 J 

-8.8128722093 J 
8.8128722093 J 
-3.0C17839895 J 
3.0017839895 J 



IMAGINARY 
0.0 J 

0.0 J 

-9.0382447685 J 
9.0382447685 J 



GAIN CONSTANT = -0.3623618484 



SYSTEM STATE VARIABLES (PHASE FORM) 



VECTOR 


1 


20863. 16713 


2 


6906. 430945 


3 


1994. 127061 


4 


124.3802348 


5 


22.38117861 


VECTOR 


5 


-0. 3623618484 


VECTOR 


1 


-283455. 3928 


2 


-27855.70425 


3 


-2180.940891 


4 


-442. 1711881 


5 


1.000000000 
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IDENTIFICATION OF UNKNOWN SYSTEM 



INITIAL 

1 

2 

3 

4 

5 



CONDITION (G) VECTOR 
16,99753824 
388.4269451 
1684. 166654 
307C4. 56626 
95381. 34774 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



3,00000 SEC 
0,50000 MSC 
120 

3 DECIMAL PLACES 
PIECEWISE CONSTANT 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 

POLES REAL 

1 -20.1316142373 

2 -1.9997989237 

3 -1.9997989237 

ZEROES REAL 

1 -8.0335677123 



IMAGINARY 
0.0 J 

-3.0011201212 J 
3.0011201212 J 

IMAGINARY 
0.0 J 



GAIN CONSTANT = 160.4525451660 



SYSTEM STATE VARIABLES (PHASE FORM) 

A VECTOR 
1 
2 
3 

B VECTOR 

3 160.4525452 

C VECTOR 

1 8.033567712 

2 l.OOOOCOOOO 

INITIAL CONDITION (G) VECTOR 

1 16.99920347 

2 418.2532176 

3 1168.376804 



261.8301183 
93, 52427869 
24. 13121208 
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EXAMPLE 



6 



SYSTEM TO BE IDENTIFIED 
SYSTEM TRANSFER FUNCTION 



POLES 


REAL IMAGINARY 


1 


-3.0000000000 0.0 J 


2 


-45.000000CG00 0.0 J 


ZEROES 


REAL IMAGINARY 


1 


-15.0000000000 0.0 J 


GAIN CONSTANT 


= 10.0000000000 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 




1 


135.0000000000 


2 


48.0000000000 


B VECTOR 




2 


10,0000000000 


C VECTOR 




1 


15.0000000000 


2 


l.OOOCOOCOOO 

% 



INITIAL STATE VECTOR 

1 1.0000000000 

2 0.0 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER 


FUNCTION 




POLES 


REAL 


IMAGINARY 


1 


-44.9825756448 


0.0 J 


2 


0. 8248831729 


-16.880278213C J 


3 


0. 8248831729 


16.8802782130 J 


4 


-2.9985948695 


0.0 J 


5 


-4.9097021934 


0.0 J 


ZEROES 


REAL 


IMAGINARY 


1 


0.8362974237 


-16.8739835601 J 


2 


0.8362974237 


16.8739835601 J 


3 


-15.0390073803 


0.0 J 


4 


-4.901716C580 


0.0 J 


GAIN CONSTANT 


= 9.9892787933 





SYSTEM STATE 


VARIABLES (PHASE FORM) 


A VECTOR 


1 


189152.5939 


2 


104719.1699 


3 


15157,98857 


4 


568. 8244217 


5 


51.24110636 


B VECTOR 


5 


9.989278793 * 


C VECTOR 


1 


21041.07999 


2 


5568.396359 


3 


325.7949073 


4 


18.26812859 


5 


l.OCOOCOOOC 
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IDENTIFICATION OF UNKNOWN SYSTEM 



INITIAL 

1 

2 

3 

4 

5 



CONDITION (G) VECTOR 
14.99958219 
633.6563568 
6068.388958 
183314.0231 
819995.9356 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



= 1.33333 SEC 
= 0.22222 MSC 
= 120 

= 4 DECIMAL PLACES 
= PIECEWISE CONSTANT 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 

POLES REAL 

1 -2.9998001615 

2 -A4. 9468936931 



ZEROES 

1 



REAL 

-14.9948012787 



IMAGINARY 
0.0 J 

0.0 J 

IMAGINARY 
0.0 J 



GAIN CONSTANT = 9.9912204742 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 134.8316990 

2 47.94669385 

B VECTOR 

2 9.991220474 

C VECTOR 

1 14.99480128 

2 1.000000000 



INITIAL CONDITION (G) VECTOR 

1 14.99876932 

2 584.3008098 
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EXAMPLE 
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SYSTEM TO BE IDENTIFIED 
SYSTEM TRANSFER FUNCTION 



POLES 


REAL IMAGINARY 


1 


-1.0000000000 -l.OCCCOOCOOO J 


2 


-1.0000000000 1.0000000000 J 


3 


-50.0COOOCOOOO 0.0 J 


ZEROES 


REAL IMAGINARY 


1 


-5.0CCOOOOOCO 0.0 J 


2 


-20.0000000000 0.0 J 


GAIN CONSTANT 


= 5.000000000C 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 




1 


100.0000000000 


2 


102.00000COOOO 


3 


52.0000000000’ 


B VECTOR 




3 


5.O00OOCOOCO 


C VECTOR 




1 


lOO.OOOCOCOOOC^ • 


2 


25.0000000000 


3 


l.OOCOOOOOOO 



INITIAL STATE VECTOR 

1 2.0000000000 



2 


l.OOOOOOOOCO 


3 


o 

• 

o 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 


1 


-49.9911920718 


o 

• 

o 


2 


1.1558667717 


-23.4137133914 


3 


1. 1558667717 


23.4137133914 


4 


-C.99972239C9 


-0.9974653718 


5 


-0.9997223909 


0. 9974653718 


ZEROES 


REAL 


IMAGINARY 


1 


1.2179838023 


-23.3753514222 


2 


1.2179838023 


23.3753514222 


3 


-19.8075068598 


o 

• 

o 


4 


-5. 1152815779 


c.o 


GAIN CONSTANT 


= 4.9560489655 





SYSTEM STATE 


VARIABLES (PHASE FORM) 


A 


VECTOR 






1 


54789.78240 




2 


55794. 37204 




3 


28434. 85333 




4 


531.2985138 




5 


49.67890331 


B 


VECTOR 






5 


4.956048965 ‘ 


C 


VECTOR 






1 


55512. 80354 




2 


13408. 14537 




3 


588.5004083 




4 


22.48682083 




5 


1 .CCOOOOOOO 



IDENTIFICATION OF UNKNOWN SYSTEM 



INITIAL 

1 

2 

3 

4 

5 



CONDITION (G) VECTOR 
225,0001273 
10975. 84397 
117558.8284 
6270060.545 
11259586.57 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



1.20000 SEC 
0.2C00C MSC 
120 

4 DECIMAL PLACES 
PIECEWISE CONSTANT 
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IDENTIFICATION OF UNKNOWN SYSTEM 
SYSTEM TRANSFER FUNCTION 

J 
J 
J 

J 
J. 

GAIN CONSTANT = 4.9685173035 



POLES 


REAL 


IMAGINARY 


1 


-49.9862042435 


0.0 


2 


-0.9995303000 


-1.0QC7669686 


3 


-0.9995303C00 


1.0007669686 


ZEROES 


REAL 


IMAGINARY 


1 


-4.9622611725 


o 

• 

o 


2 


-20.2677062457 


o 

• 

o 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 100.0021676 

2 101.9260468 

3 51.98526484 



B VECTOR 

3 4.968517303 



C VECTOR 

1 100.5736518 

2 25.22996742 

3 l.OOOCOOOOO 



INITIAL CONDITION (G) VECTOR 

1 225.0005829 

2 11494.69251 

3 20485.03332 
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SYSTEM TO BE IDENTIFIED 



SYSTEM TRANSFER 

POLES 

1 

2 

3 

4 



NOTION 

REAL 

-3.0000000000 

- 8.0000000000 

-40.0000000000 

-40.0000000000 



IMAGINARY 
0.0 J 

0.0 J 

-5.0000000000 J 
5.00C00C0000 J 



ZEROES 

1 



REAL 

- 20.0000000000 



IMAGINARY 
0.0 J 



GAIN CONSTANT = 160.0000000000 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 
1 
2 

3 

4 



39000.0030000000 
19795. OOOOOCOOOO 
2529.00OOOCOOO0 
91.0000000000 



B VECTOR 
4 



160.0000000000 



C VECTOR 
1 
2 



20.0000000000 

l.OOOOOOOOoO 



INITIAL STATE VECTOR 

1 3.0000COOOOO 

2 2.0000000000 

3 1.0000000000 

4 0.0 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL 


IMAGINARY 


1 


-40.0576848888 


-4. 8648939727 


2 


-40.0576848888 


4. 8648939727 


3 


-8.0002609098 


o 

• 

o 


4 


-1.4530133153 


o 

• 

o 


5 


-2.9998483817 


o 

• 

o 


ZEROES 


REAL 


IMAGINARY 


1 


-1.4492282436 


o 

• 

o 


2 


-20.3685768935 


0.0 


GAIN CONSTANT 


= 157.7509613037 





SYSTEM STATE 


VARIABLES (PHASE FORM) 


A 


VECTOR 






1 


56781.06779 




2 


67897.28681 




3 


23515.35114 




4 


2665. 954710 




5 


92. 56849238 


B 


VECTOR 






5 


1 57.7509613 


C 


VECTOR 


k ■ 




1 


29.51871692 




2 


21.81780514 




3 


1.000000000 
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IDENTIFICATION OF UNKNOWN SYSTEM 



INITIAL 

1 

2 

3 

4 

5 



CONDITION (G) VECTOR 
62,00108720 
5780,206688 
169108,2337 
1409764,537 
1708815,026 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



1,48842 SEC 
0,24807 MSC 
120 

5 DECIMAL PLACES 
PIECEWISE CONSTANT 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER 

POLES 

1 

2 

3 

4 



FUNCTION 

REAL 

-40.0840542053 

-40.0840542053 

-3.0000451598 

-7.99971C0927 



ZEROES 

1 



REAL 

-19.9584019372 



IMAGINARY 
-4.7409985601 J 
4.7409985601 J 
0.0 J 

0.0 J 

IMAGINARY 
0.0 J. 



GAIN CONSTANT = 160.72850C3662 



SYSTEM STATE VARIABLES (PHASE FORM) 

A VECTOR 
1 
2 

3 

4 

B VECTOR 
4 

C VECTOR 

1 19.95840194 

2 1.000000000 

INITIAL CONDITION (G) VECTOR 

1 62.00114277 

2 5693.349573 

3 160935.4623 

4 1176687.962 



39100.17487 
19844. 88825 
2535.037532 
91.16786366 



16C. 7285004 
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SYSTEM TO BE IDENTIFIED 
SYSTEM TRANSFER FUNCTION 



POLES 


REAL IMAGINARY 


1 


-1.0000000000 0.0 J 


2 


-5.0000000000 0.0 J 


ZEROES 


REAL IMAGINARY 


GAIN CONSTANT = 


lO.OOOOOOOOOC 



SYSTEM STATE VARIABLES (PHASE FORM) 
A VECTOR 



1 


5.0000000000 


2 


6.0000000000 


B VECTOR 





2 lO.OCOOOOOOCO 



C VECTOR 




1 


1.0000000000 



INITIAL STATE VECTOR 



1 


l.OOOOOCOOOO 


2 


o 

• 

o 
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IDENTIFICATION OF UNKNOWN SYSTEM 
SYSTEM TRANSFER FUNCTION 

J 
J 
J 
J 
J 



J 
J 

GAIN CONSTANT = 28777.0156250000 



POLES 


REAL 


IMAGINARY 


1 


-5.C110029817 


o 

• 

o 


2 


-2874. 4029427054 


o 

• 

o 


3 


C.C094513173 


-2.4186034399 


4 


0.0094513173 


2.4186034399 


5 


-0.9999781363 


o 

• 

o 


ZEROES 


REAL 


IMAGINARY 


1 


0.0096691675 


-2.4198473522 


2 


0.0096691675 


2.4198473522 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 84255.60063 

2 100828.6134 

3 30^)26.28200 

4 17234.39503 

5 2880.395021 



B VECTOR 

5 28777.01563 



C VECTOR 

1 5.855754701 

2 -0. 1933833505D-01 

3 1.000000000 
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IDENTIFICATION OF UNKNOWN SYSTEM 



INITIAL 

1 

2 

3 

4 

5 



CONDITION (G) VECTOR 
1.450317645 
2876.228032 
17252. 85740 
16462.65716 
101232.6568 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



= 12.00000 SEC 
= 2.000C0 MSC 
= 120 

= 5 DECIMAL PLACES 
= PIECEWISE CONSTANT 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 



POLES 


REAL IMAGINARY 


1 


-0.9999961462 0.0 J 


2 


-5.0000869786 0.0 J 


ZEROES 


REAL IMAGINARY 


GAIN CONSTANT 


= 10.0001306534 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 




1 


5.000067709 



2 6.000083125 



B VECTOR 




2 


10.00013065 


C VECTOR 




1 


1.000000000 



INITIAL CONDITION (G) VECTOR 

1 1.00002578A 

2 6.000036891 

% 
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EXAMPLE 10 



SYSTEM TO BE IDENTIFIED 



SYSTEM TRANSFER FUNCTION 

POLES REAL 

1 -5.00000OC0OO 

2 -5.0000000000 

3 -100.0000000000 



IMAGINARY 
-5.0000000000 J 
5.000000C000 J 
0.0 J 



ZEROES 

1 



REAL 

- 20.0000000000 



IMAGINARY 
0.0 J 



GAIN CONSTANT = 100. OCOOCOOOOO 



SYSTEM STATE VARIABLES (PHASE FORM) 
A VECTOR 

1 5000.00000C0000 

2 1050.0000000000 

3 110.0000000000 



B VECTOR 
3 



100.0000000000 



C VECTOR 

1 20.0000000000 
2 1 . 0000000000 » 



INITIAL STATE VECTOR 

1 2.0000000000 

2 1.0000000000 

3 0.0 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER FUNCTION 

POLES REAL 

1 3533.8075084749 

2 -25.9614103219 

3 -100.0395044818 

4 -4.9998977161 

5 -4.9998977161 



ZEROES REAL 

1 -19.5630088534 

2 -26.9189879242 



IMAGINARY 
0.0 J 

0.0 J 

0.0 J 

-5.0000893933 J 
5.000C893933 J 



IMAGINARY 
C.O J 

0.0 J 



GAIN CONSTANT =-348534.3750000000 



SYSTEM STATE VARIABLES (PHASE FORM) 



A VECTOR 

1 -458893162.9 

2 -11391G225.2 

3 -13774844.37 

4 -476693.1811 

5 -3397.806798 

B VECTOR 

5 -348534.3750 



C VECTOR 

1 526.6163991 

2 46.48199678 

3 1.000000000 
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IDENTIFICATION OF UNKNOWN SYSTEM 



INITIAL 

1 

2 

3 

4 

5 



CONDITION (G) VECTOR 
40, 97649863 
-139289.1252 
-19623611.76 
-535751377.4 
-3138782043. 



PROGRAM PARAMETERS 
RECORD LEGNTH 
SAMPLING PERIOD 
EQUATIONS FORMED 
PRECISION OF DATA 
INPUT FUNCTION 



0,600GC SEC 
C.IOCCO MSC 
120 

5 DECIMAL PLACES 
PIECEWISE CONSTANT 
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IDENTIFICATION OF UNKNOWN SYSTEM 



SYSTEM TRANSFER 

POLES 

1 

2 

3 

ZEROES 

1 



FUNCTION 

REAL 

-99,9959669612 

-4.9999968528 

-4.9999968528 

REAL 

-19,9960841378 



GAIN CONSTANT = 100.0136718750 



IMAGINARY 
0,0 J 

-5,0000354476 J 
5.0000354476 J 

IMAGINARY 
0.0 J 



SYSTEM STATE VARIABLES (PHASE FORM) 
A VECTOR 

1 4999.830647 

2 1049,959363 

3 109.9959607 

B VECTOR 

3 100.0136719 



C VECTOR 

1 19,99608414 

2 l.COOOOOOOO 



INITIAL CONDITION (G) VECTOR 

1 40.99999105 

2 4529.827485 

3 34198. 63350 



7 4 



IV. CONCLUSIONS 



The method of multiple integrations is a practical 
3-ud flexible method for identifying lumped parameter, 
time invariant systems. Complete and accurate 
identifications can be made on the basis of arbitrary in- 
put-output records taken over a short time interval and 
accurate to only three or four significant figures. The 
computational requirements of the method are not exces- 
The procedure can be implemented by relying en- 
tirely on subroutines which are available in most computer 
center libraries. 

At present the technique is limited to comparatively 
low order systems. When the input-output records are good 
to three or four significant digits the method will be 
capable of identifying systems up to about fifth order. 
This limitation is due primarily to the algorithm used to 
solve the overdetermined set of linear equations. As 
better algorithms become available it will be possible to 
identify higher order systems. 

The accuracy of an identification is usually compar- 
able to the number of significant digits in the input- 
output data. Accuracy depends to a lesser extent on the 
sampling rate, the nature of the input function driving 
the system, and the order of the system. 

There are several areas where additional research 
might prove fruitful. Since the input-output records must 



7S 




f 




be in sampled data form for the computer it might be prof- 
itable to reformulate the identification procedure from a 
sampled system standpoint. Standard techniques could be 
used to convert the sampled system representation obtained 
from the identification to a continuous system representa- 
tion. The ability to identify sampled systems would be a 
worthwhile extension of this method. 

The method of multiple integrations may prove very 
useful as a tool for approximating high order systems with 
low order systems. Research could be done to determine 
the quality of the approximations obtained using this 
method. 
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DIGITAL COMPUTER PROGRAMS 



MAIN PROGRAM - LINEAR SYSTEM IDENTIFICATION 



PURPOSE 

TO IDENTIFY LINEAR TIME INVARIANT 
THE BASIS OF INPUT-OUTPUT RECORDS 



SYSTEMS ON 



DESCRIPTION 

INPUT 

NP 

NZ 

KPTMAX - 

IPTS 

T 

R 

C 



OF PARAMETERS 



ESTIMATED NUMBER OF POLES 
ESTIMATED NUMBER OF ZEROES 
NUMBER OF DATA POINTS 
DATA POINTS INTEGRATED PER 
TIME 

INPUT AMPLITUDE AT TIME T 
OUTPUT AMPLITUDE AT TIME T 



LINEAR EQ. 



SIMULATED SYSTEMS 



REMARKS 

(1) OUTPUT WILL CONSIST OF A TRANSFER FUNCTION 
AND STATE VARIABLE REPRESENTATION OF SYSTEM 

(2) PROGRAM IS PRESENTLY CONFIGURED TO IDENTIFY 
SYSTEMS SIMULATED BY SUBROUTINE SYSTEM. IF 
IT IS DESIRED TO IDENTIFY A PHYSICAL SYSTEM 
REPLACE SUBROUTINE SYSTEM CALLS WITH APPROP- 
RIATE READ STATEMENTS. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
WHEN IDENTIFY/ING PHYSICAL SYSTEMS 

(1) DLLSQ 

(2) RTPLSB 
WHEN IDENTIFYING 

(1) DLLSQ 

(2) RTPLSB 

(3) SYSTEM 
(A) EXPAND 

(5) SUMM 

(6) DIFF 

(7) PROD 

(8) GAUSS3 

(9) ARRAY 

(10) MINV 

(11) ROUND 

METHOD 

MULTIPLE INTEGRALS OF THE INPUT AND OUTPUT DATA 
ARE USED TO FORMULATE A SET OF OVERDETERMINED 
LINEAR EQUATIONS. THESE EQUATIONS ARE THEN SOLVED 
FOR THE UNKNOWN MODEL PARAMETERS USING THE METHOD 
OF LEAST SQUARES. 



5 



REALM'S T0,TN(11),R0(11),RN( 
REAL*8 A(26C0) ,8(200) , X(26) 
REAL*8 AP(IC), AZ(10),PRA(9) 
INTEGER IPIV(26) 

KPTMAX=5002 

IPTS=5C 

NP=4 

NZ=NP-1 

CONTINUE 

MEQS=KPTMAX/IPTS 
T0=-1.00 
EPS=10.0*^(-35 ) 

K=1 

M=0 

NP1=NP+1 

NP2=NP+2 



11 ) , COdl ) ,CN( 
,AUX(?2) ,CONV( 
,PIA(9),ZRA(8) 



11) ,DT2 
9) 

, ZIA(8) 
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NZ1=NZ+1 

N=NP+NZ+1 

NI=N+NP 

SET CUMULATIVE INTEGRAL VALUES TO ZERO 

DO 10 1=2, NP2 
R0( I )=0.06 
10 C0(I)=0.0C 

READ IN INITIAL DATA POINT (T,R,C) 

CALL SYSTEMITO, R0(1 ) ,C0( 1 ) ) 

T0FF=T0 
15 K=K+1 

READ IN NEW DATA POINT (T,R,C) 

CALL SYSTEMITNd) ,RN(1) ,CN(1) ) 

UPDATE MULTIPLE INTEGRATIONS 

DT2=(TN( 1)~T0)’:'0.5 
DO 20 INT=1,NP1 

RN( INT+1 )=(R0( I NT)+RN( INT) )*DT2+R0( INT+1 ) 

20 CN ( INT+1 ) = (C0( INT)+CN( INT) ) =«'DT2+C0( I NT+ 1 ) 

FORM A LINEAR EQUATION 

IF (K.NE. (K/IPTS )#IPTS) GO TO 35 
M=M+1 

B(M)=CN( 2) 

TN(2)=(TN(1 )-T0FF) 

DO 25 1=1, NP 
I A=( NP-I )+-MEOS+M 
IC=(N+I-1 )*MEQS+M 

TN( I+2)=TN( I+l )*TN( 2) /FLOAT! I+l) 

A(IA)=-CN( 1+2) 

25 A( IC)=TN( I+l ) 

DO 30 1=1, NZl 
IA=(NP+I-1 )*MEQS+M 
IRN=NP+3-I 
30 A( IA)=RN( IRN) 

RESET OLD VALUES 

35 T0=TN( 1 ) 

DO 40 1=1, NP2 
R0( I )=RN( I ) 

40 C0( I )=CN( I ) 

IF (M.LT.MEQS) GO TO 15 

SOLVE FOR PARAMETERS BY METHOD OF LEAST SQUARES 

CALL DLLSQ( A,B,M,NI ,1, X, IPiViEPS, IER,AUX) 

CALCULATE POLES 

AP(1)=1.C0 
DO 45 1=1, NP 
J=NP+2-I 
45 AP ( J ) =X ( I ) 

CALL RTPLSB(NP,AP, PRA,PIA,CONV, lERPZ) 

CALCULATE ZEROES 

GAINI=X(N) 

DO 50 I=NP1,N 
50 X( I )=X( I )/X!N) 

DO 55 1=1, NZl 
J=N+1-I 
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55 /iZ(I) = X(J) 

IF (NZ.EQ.O) GO TO 60 

ZRA, ZIAtCONV, lERPZ) 

60 CONTINUE 
OUTPUT 

WRITE(6,914) 

WRITE(6,915) K 
WRITE(6,916) M 
WRITE(6,917) IPTS 
WRITE(6,918) EPS 
WRITE(6,919) AUX(l) 

WRITE(6,920) lER 
WRITE(6,921) ( IPIV( I) »I=1,NI ) 

WRITE(6,900) 

WRITE(6t901 ) 

WRITE(6,902) 

DO 65 1=1, NP 

WRITE(6,903) I , PRA ( I ) , P I A ( I ) 

65 CONTINUE 

WRITE(6,904) 

IF (NZ.EQ.C) GO TO 75 

00 70 1=1, NZ 

WRITE(6,903) I , ZRA ( I ) , Z I A ( I ) 

70 CONTINUE 
75 CONTINUE 

WRITE(6,905) GAINI 
WRITE(6,906) 

WRITE(6,9C7) 

DO 80 1=1, NP 

1 I = NP+2-I 

WRITE(6,91C) I,AP(II) 

80 CONTINUE 

WRITE(6,908) 

WR!TE(6-910) NP, GAINI 
WRITE(6,909) 

DO 85 1=1, NZl 
I I=NZ+2-I 

WRITE(6,910) I,AZ(II) 

85 CONTINUE 

WRITE(6,913) 

DO 90 1=1, NP 
I I = N + I 

WRITE(6,910) I,Xni) 

90 CONTINUE 
GO TO 5 

900 F0RKAT(1H1,///,25X, 'IDENTIFICATION OF UNKNOWN SYSTEM') 

901 FORMAT!///, 12X, 'SYSTEM TRANSFER FUNCTION') 

902 FORMAT (//, 15X, ' POLES' , 12X, 'REAL ' , 13X , • 1 M AG I NARY • , / ) 

903 F0RMAT(17X,I2,7X,G15.8,6X,G15.8,1X, 'J',/) 

904 FORMAT! //, 1 5X, ' ZEROES ', 1 1 X ,' RE AL ', 1 3X , ' IMAGINARY* ,/ ) 

905 FORMAT!//, 15X, 'GAIN CONSTANT =',G15.8,/) 

906 FORMAT!////, 12X, 'SYSTEM STATE VARIABLES IPHASE FORM)') 

907 FORMAT! //, 15X, ' A VECTOR',/) 

9C8 FORMAT!//, 15X, ' B VECTOR',/) 

9C9 FORMAT!//, 15X, 'C VECTOR',/) 

910 F0RMATI17X, I2,7X,G15.8,/) 

911 F0RMAT!1H , 9X , ' I ER ' , I 5 , 5X , ' E PS ' , E 1 5 . 8 , 5X , ' AUX ' , 
1E16.8,/,9X, 111 8) 

912 F0RMAT!//,5X,4!E20.9,5X) ,//, 5X,4!E20.9,5X) ,//) 

913 FORMAT!//, 15X, ' INITIAL STATE VECTOR',/) 

914 FORMAT !//////, 12X, ' PROGRAM PARAMETERS’,/) 



915 F0RMAT!15X, 

916 F0RMATI15X, 

917 FORMAT! 15X, 

918 F0RMAT!15X, 

919 FORMAT! 15X, 

920 FORMAT! 15X, 

921 FGRMAT!15X, 
STOP 

END 



NUMBER OF DATA POINTS =',I6,/) 
NUMBER OF EQUATIONS =',I6,/) 
DATA POINTS PER EQUATION =',I6,/) 



EPS 



=' , E16.8 ,/ ) 



RMS ERROR =' ,E16.8,/) 
lER =',I4,/) 

IPIVm =' ,15! I2,1X) ,/) 
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SUBROUTINE SYSTEM 



PURPOSE 

GENERATES INPUT-OUTPUT RECORDS FOR IDENTIFICA- 
TION PROGRAM BY SIMULATING A SYSTEM DESCRIBED 
BY A TRANSFER FUNCTION THAT IS READ IN 

USAGE 

CALL SYSTEM(T,R,C) 

DESCRIPTION OF PARAMETERS 
INPUT 

NP - NUMBER OF POLES 

P(I) - THE VECTOR OF POLES 
NZ - NUMBER OF ZEROES 

Z(I) - THE VECTOR OF ZEROES 
GAIN - THE GAIN CONSTANT 
OUTPUT 

T - TIME 

R - INPUT AMPLITUDE AT TIME T 

C - OUTPUT AMPLITUDE AT TIME T 

REMARKS 

(1) INPUT FUNCTION MAY BE CHANGED BY 
CHANGING ONE CARD IN PROGRAM 

(2) R AND C ARE ROUNDED TO NA DIGITS 

(3) ALL FLOATING POINT VARIABLES ARE DECLARED 
DOUBLE PRECISION (REAL-'B). 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

(1) SUMM 

(2) DIFF 
(3» PROD 

(4) EXPAND 

(5) GAUSS3 

(6) ARRAY 

- (7) MINV 

(8) ROUND 

METHOD 

TRANSFER FUNCTION IS CONVERTED TO STATE 
VARIABLE REPRESENTATION AND INTEGRATED USING 
TRAPEZOIDAL INTEGRATION 



SUBROUTINE SYSTEMIT , R,C) 



REAL*8 AA( 10) ,CC( 10) ,A(9,9> ,B(9t 1 ) , XXX( 9,9) 
REALMS PHI (9,9 ) , DEL (9, 1 ) , XX(9, 1 ) ,UU( 1 , 1 ) 
REAL’!' 8 T,R,C,U,DT,AI(9,9),ZZZ(9,9) 
COMPLEX’i'16 P(9) ,Z( 8) 



IF (T.GE.0.00) GO TO 55 

INPUT POLES, ZEROES, AND GAIN CONSTANT 



5 

10 



^EAD(5,899,END=999) NP 



^EAD(5,898) 

:ALL EXPAND 
DO 5 1 = 1, NP 

:c( I ) =0.00 

DONTINUE 
^EAD(5, 899) 

IF (NZ.EQ.O) 
^EAD(5,898) 

:ALL EXPAND 
:ONTINUE 
^EAD(5,898) GAIN 



(P( I ) , 1=1 ,NP) 
(NP, P, AA) 



NZ 

GO TO ID 
(Z( I ) , 1=1, 
(NZ,Z,CC) 



NZ) 
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SET NA 
NA=3 

FIND SHORTEST TIME CONSTANT AND CALCULATE DT 

DT=0.00 
DO 15 1=1, NP 

IF { DT.LT.CDABS( P( I ) ) ) DT=CDABS ( P ( I ) ) 

15 CONTINUE 

IF (NZ« EQ.O) GO TO 25 
DO 20 1=1, NZ 

IF (DT.LT.CDABS(Z(I) )) DT=CDABS ( Z (I ) ) 

20 CONTINUE 
25 DT=1.0/(DT*100.0) 

IF (DT. EQ.O. 00) DT=0.00C01 

FORM A, B, AND C MATRICES 

DO 35 1=1, NP 
DO 30 J=1,NP 
AI( I , J) =0.00 
A( I, J)=0.00 
30 CONTINUE 
B( I , 1)=0.00 
35 CONTINUE 

DO 40 1=2, NP 
All I , I ) = 1.00 
A( I-l, ! ) = DT/2.00 
A( NP, I )=-AA( I ) *DT/2.00 
40 CONTINUE 

AI(1,1)=1.00 
A(NP,1 )=-AA(l )=!'DT/2.00 
B(NP, 1) =GAI NvOT/2.00 
CC(NZ+1 )=i.on 

CALCULATE PHI MATRIX 

CALL DIFF(A1,A,ZZZ,NP,MP) 

CALL GAUSS3(NP , EPSS ,Z Z Z , XXX , KER , 9 ) 

CALL SUMM(AI,A,ZZZ,NP,NP) 

CALL PROD(XXX,ZZZ,PHI ,NP,NP,NP) 

CALCULATE DEL MATRIX 

CALL PR0D(XXX, B,0EL,NP,NP, 1) 

DEFINE INITIAL CONDITIONS 

K = 1 

DO 45 I=1,NP 

XX( I , 1 )=FL0AT{ NP-I ) 

45 CONTINUE 

UU( 1 , 1 ) =0.00 
T=0.00 
R=0.00 
C=0.00 

DO 50 1=1, NP 
C=CC( I )*XX( I ,1 )+C 
50 CONTINUE 
GO TO 65 

CALCULATE NEW DATA POINT (T,R,C) 

55 CONTINUE 

T=DT*FLOAT(K) 

INPUT FUNCTION 

U= +0.1*FL0AT(K/53)+FL0AT(K/4C3)+FL0AT(K/6C3) 
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UU(1,1 )=UU(1,1 )+U 

CALL PR0D(PHI,XX,XXX,NP,NP,1) 

CALL PRODCDEL, UU,ZZZ,NP, 1 ,1) 

C — 0 • 00 

DG 60 1=1, NP 

XX( I , 1 )=XXX( I, 1 )+ZZZ( I ,1) 

C=CC( I )*XX( 1,1 )+C 
60 CGNTINUE 
R=U 

ROUND OFF R AND C 

CALL R0UND(R,NA) 

CALL R0UND(C,NA) 

UU(1,1)=U 

K=K+1 

RETURN 

OUTPUT 

CONTINUE 
WRITE(6,900) 

WRITE(6,901 ) 

WRITE(6,902) 

DO 70 1=1, NP 
WRITE(6,903) I,P(I) 

CONTINUE 
WRITE(6,90A) 

IF (NZ.EQ.O) GO TO 80 
DO 75 1=1, NZ 
WRITE(6,903) I,Z(I) 

CONTINUE 
CONTINUE 

WRITE(6,905) GAIN 
WRITE(6,906) 

WRITE(6,907) 

DO 85 1=1, NP 
WRITE(6,910) I,AA(I) 

CONTINUE 
WRITE(6,908) 

WRITE(6,910) NP,GAIN 
WRITE (6,909) 

NZ1=NZ+1 
DO 90 1=1, NZl 
WRITE(6,910) I,CC(I) 

CONTINUE 
WRITE(6,913 ) 

DO 95 1=1, NP 
WRITE(6,910) I,XX(I,1) 

CONTINUE 

FGRM.AT( 2F10.5) 

FORMAT( n ) 

FORMAT( 1H1,////,25X, ’SYSTEM TO BE IDENTIFIED') 

FORMAT! ///, 12X, ' SYSTEM TRANSFER FUNCTION') 

FORMAT ( //, 15X, ' POLES' , 12X, • REAL • , 13X , • IMAGINARY* , / ) 
FORMAT! 17X,I2,7X,F]4.7,6X,F1A.7,1X,'J',/) 

FORMAT! //,15X, ' ZEROES' ,1 IX, 'REAL' ,13X,' IMAGINARY' ,/) 
FORMAT! //, 15X, ' GAI N CONSTANT =',F14.7,/) 

FORMAT!////, 12X, 'SYSTEM STATE VARIABLES (PHASE FORM)') 
FORMAT! //, 15X, ' A VECTOR',/) 

FORMAT! //,15X, ' B VECTOR',/) 

FORMAT! //, 15X, ' C VECTOR',/) 

FORMAT! 17X, 12, 7X,F14.7,/) 

FORMAT!//, 15X, ' INITIAL STATE VECTOR’,/) 

RETURN 
CONTINUE 
STOP 
END 



65 



70 



75 

80 



85 



90 



95 

898 
8 99 

900 

901 

902 

903 

904 

905 

906 

907 

908 

909 

910 
913 



999 
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c 
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SUBROUTINE DLLSQ 



PURPOSE 

TO SOLVE LINEAR LEAST SQUARES PROBLEMS, I.E. TO 
MINIMIZE THE EUCLIDEAN NORM OF B-A<=X, WHERE A IS 
A M BY N MATRIX WITH M NOT LESS THAN N. IN THE 
SPECIAL CASE M=N SYSTEMS OF LINEAR EQUATIONS MAY 
BE SOLVED. 

USAGE 

CALL DLLSQ(A,B,M,N,L,X,I PIV,EPS,IER,AUX) 



DESCRIPTION 

A 



M 

N 



L 

X 

IPIV 



EPS 



lER 

AUX 



OF PARAMETERS 

DOUBLE PRECISION M BY N MATRIX 
(DESTROYED). 

DOUBLE PRECISION M BY L RIGHT HAND SIDE 
MATRIX (DESTROYED). 

ROW NUMBER OF MATRICES A AND B 

COLUMN NUMBER OF MATRIX A, ROW NUMBER OF 

MATRIX X 

COLUMN NUMBER OF MATRICES B AND X 
DOUBLE PRECISION N BY L SOLUTION MATRIX 
INTEGER OUTPUT VECTOR OF DIMENSION N 
WHICH CONTAINS INFORMATION GN COLUMN 
INTERCHANGES IN MATRIX A. 

SINGLE PRECISION INPUT PARAMETER WHICH 
SPECIFIES A RELATIVE TOLERANCE FOR 
DETERMINATION OF RANK OF A. 

A RESULTING ERROR PARAMETER 
A double precision AUXILIARY STORAGE 
ARRAY OF DIMENSION MAX(2^^M,L). ON RETURN 
FIRST L LOCATIONS OF AUX CONTAIN THE 
RESULTING LEAST SQUARES. 



REMARKS 

(1) NO ACTION BESIDES ERROR MESSAGE IER=-2 IN 
CASE M LESS THAN N. 

(2) NO ACTION BESIDES ERROR MESSAGE IER=-1 IN 
CASE OF A ZERO MATRIX A. 

(3) IF RANK K OF MATRIX A IS FOUND TO BE LESS 
THAN N BUT GREATER THAN C , THE PROCEDURE 
RETURNS WITH ERROR CODE IER=K INTO CALLING 
PROGRAM. THE LAST N-K ELEMENTS OF VECTOR IPIV 
DENOTE THE USELESS COLUMNS IN MATRIX A. 

(4) IF THE PROCEDURE WAS SUCCESSFUL, ERROR 
PARAMETER lER IS SET TO ZERO. 



SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 



METHOD 

HOUSEHOLDER 
MATRIX A TO 
APPLIED THE 
APPROXI MATE 
BY BACK SUBSTITUTION 
GOLUB, G. , NUMERICAL 



TRANSFORMATIONS ARE USED TO TRANSFORM 
UPPER TRIANGULAR FORM. AFTER HAVING 
SAME TRANSFORMATIONS TO MATRIX B, AN 
SOLUTION OF THE PROBLEM IS COMPUTED 
FOR REFERANCE, SEE 
METHODS FOR SOLVING LINEAR 



LEAST SQUARES PROBLEMS, NUMERISCHE 
VOL. M, ISS.3 (1965), PP. 206-216. 



MATHEMATIK. 



SUBROUTINE DLLSQ ( A , B, M , N, L , X , IPIV, EPS, lER, AUX) 



DIMENSION A(1),B(1),X(1),IPIV(1),AUX(1) 

DOUBLE PRECISION A , B , X , AUX , P I V, H , S I G , BETA , TOL 



ERROR TEST 
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IF(M-N) 30, 1, 1 

GENERATICN OF INITIAL VECTOR S(K) (K=l,2,..., 
IN STORAGE LOCATIONS AUX(K) ( K= 1 , 2 , . . . , N ) 



N) 



1 PIV=0.D0 
IEND=0 

DC 4 K=1,N 
IPIV(K) =K 
H=0.D0 
IST=IEND+1 
I END=IEND+M 
DO 2 I=IST,IEND 

2 H=H+A( I )*A( I ) 

AUX(K)=H 

IF(H-PIV)4,4,3 

3 PIV=H 
KPIV=K 

4 CONTINUE 

ERROR TEST 
IF(PIV)31,3L,5 

DEFINE TOLERANCE FOR CHECKING RANK OF A 

5 SIG=DSQRT( PIV) 

T0L = SIG=!=ABS(EPS) 



DECOMPOSITI CN LOOP 

LP=L=^^M 

IST=-M 

DO 21 K=1,N 

IST=IST+M+1 

I END=I ST+M-K 

I=^KPIV-K 

IF( I )8,8,6 

INTERCHANGE K-TH COLUMN OF A WITH KPIV-TH IN CASE 
KPIV.GT.K. 

6 H=AUX(K) 

AUX(K)=AUX( KPI V) 

AUX( KPIV )=H 
ID=I*M 

DO 7 I=IST,IEND 

J=I+ID 

H=A( I ) 

A(I )=A( J) 

7 A(J)=H 

COMPUTATION OF PARAMETER SIG 

8 IF(K-1 ) 11, 11,9 

9 SIG^^O.OO 

DO 10 I=IST,IEND 

10 SIG=SIG + A( I )'^A( I ) 

SIG=DSQRT( SIG) 

TEST ON SINGULARITY 
IF( SIG-T0D32, 32,11 

GENERATE CORRECT SIGN OF PARAMETER SIG 

11 H=A(IST) 

IF(H) 12,13,13 

12 SIG=-SIG 

SAVE INTERCHANGE INFORMATION 

13 IPIV(KPIV)=IPIV(K) 

IPIV(K)=KPIV 

GENERATION OF VECTOR UK IN K-TH COLUMN OF MATRIX 

A AND OF PARAMETER BETA 

BETA=H+SIG 
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A( IST)=BETA 
BETA=1.00/(SIG=i=BETA) 

J=N+K 

AUX{ J)=-SIG 
IF(K-N) 14tl9tl9 

TRANSFORMATION OF MATRIX A 

14 PIV=O.DO 
ID=C 
JST=K+1 
KPIV=JST 

DO 18 J=JSTtN 

ID=ID+M 

H=0.D0 

DO 15 I=IST,IEND 
II=I+ID 

15 H=H+A( I ) ’J'AC 1 1 ) 

H=BETA«H 

DO 16 I=IST,IEND 
I 1=1 + 10 

16 A{ I I )=A( II )-A( I )«H 

UPDATING OF ELEMENT S(J) STORED IN LOCATION AUX(J) 
I I = IST + ID 

H=AUX(J)-A( II )*A(II ) 

AUX{ J)=H 

IF(H-PI V)18,18tl7 

17 PIV=H 
KPIV=J 

18 CONTINUE 

TRANSFORMATION OF RIGHT HAND SIDE MATRIX B 

19 DO 21 J=K,LMtM 
H=0.00 
IEND=J+M-K 

I I=IST 

DO 20 I=J, I END 
H=H+A( I I )*B{ I ) 

20 11 = 11+1 
H=BETA=f=H 
I I=IST 

DO 21 I=J, I END 
B( I ) = B( I )~A( I I )*H 
21 11 = 11+1 

END OF DECOMPOSITION LOOP 



BACK SUBSTITUTION AND BACK INTERCHANGE 

IER=0 

I=N 

LN=L=«=N 

PIV=1.D0/AUX(2*N) 

DO 22 K=NtLNtN 
X(K)=PI V^B( I ) 

22 I=I+M 

IFIN-l 126,26,23 

23 JST=(N-1 )^M+N 
DO 25 J=2,N 
JST=JST-M-1 
K=N+N+1-J 
PIV=1.D0/AUX(K ) 

KST=K-N 

ID=IPIV(KST)-KST 
IST=2-J 
DO 25 K=1,L 
H=B(KST) 

IST=IST+N 

lEND=IST+J-2 

II=JST 

DO 24 I=IST,IEN0 
II=II+M 

24 H=H-A( I I )*X( I ) 
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I=IST-1 
II=I+ID 
X( I ) = X( I I ) 
X( II )=PIV«H 
25 KST=KST+M 



COMPUTATION OF LEAST SQUARES 

26 IST=N+1 
IEND=0 

DO 29 J=1»L 
IEND=IEND+M 
H=O.DO 

IF(M-N)29t29,27 

27 DO 28 I=IST,IEND 

28 H=H+B( I )*B( I ) 

IST=IST+M 

29 AUX(J)=H 
RETURN 

ERROR RETURN IN CASE M LESS THAN N 

30 IER=-2 
RETURN 

ERROR RETURN IN CASE OF ZERO-MATRIX A 

31 IER=-1 
RETURN 

ERROR RETURN IN CASE OF RANK OF MATRIX A LESS THAN N 

32 IER=K-1 
RETURN 
END 



SUBROUTINE ROUND 
PURPOSE 

TO ROUND OFF A NUMBER TO A SPECIFIED NUMBER OF 
SIGNIFICANT DIGITS 

USAGE 

CALL ROUND(AtN) 

DESCRIPTION OF PARAMETERS 

A - NUMBER TO BE ROUNDED 

N - SIGNIFICANT DIGITS TO BE RETAINED 

REMARKS 

NONE 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 



SUBROUTINE ROUND(AtN) 

REAL’!'8 A 

X=A 

IF (X.EQ.0.0) GO TO 1 
SIGN=+1.0 

IF (X.LT.O.O) SIGN=-1.0 
L=AL0G10(A3S(X( )+l.CO 
Y=X»10.0**(N-L ) 

I=Y 

Z=I 

IF ( (Y-Z).GT.0.5) 1=1+1 
Z = I 

A=SIGN#Z*10.0=«'*( L-N) 

1 CONTINUE 
RETURN 
END 
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SUBROUTINE RTPLSB 



PURPOSE 

TO FIND THE ROOTSt BOTH REAL AND COMPLEX, OF A 
POLYNOMIAL WITH REAL COEFFICIENTS USING BOTH 
BAIRSTOW AND NEWTON-RAPHSON METHODS. 

USAGE 

CALL RTPLSB(N,A,U,V,CONV, lER) 

DESCRIPTION OF PARAMETERS 
INPUT 

N - DEGREE OF POLYNOMIAL 

A - COEFFICIENT VECTOR OF POLYNOMIAL 

OUTPUT 

U - VECTOR OF REAL PARTS OF ROOTS 

V - VECTOR OF IMAGINARY PARTS OF ROOTS 

CONV - CONVERGENCE INDICATORS FOR EACH ROOT 

lER - ERROR INDICATOR 

=0 , N IS WITHIN BOUNDS 
=1 , N IS LESS THAN ONE 

REMARKS 

(1) FOR PROBLEMS WITH NONMULTIPLE ROOTS, THE 
APPROXIMATE NUMBER OF SIGNIFICANT DIGITS OF 
EACH PART OF EACH ROOT WILL APPEAR AS THE 
EXPONENT OF THE CORRESPONDING ENTRY IN CONV 

(2) ACCURACY MAY BE LESS THAN INDICATED BY CONV 
IF THERE ARE MULTIPLE ROOTS 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 



SUBROUTINE RTPLSB ( N , A , U , V , CONV , I ER ) 
IMPLICIT REAL’!''8 (A-H), REAL-8 (0-Z) 
DIMENSION A( 1 ) ,U(1 ) ,V(1 ) ,CONV( 1 ) , B( 53) 
DIMENSION C(53) ,D( 53) ,E(53) ,H(53) 

KF = 10 
L = 25 
IER=0 
NN=N 
K=0 

N1=N+1 

IZF=0 

IF(NN) 51,51,52 
52 DO 11 1=1, N1 

IF(A(I)) 15,12,15 
15 CONTINUE 
IZF=1 
GO TO 11 

12 CONTINUE 

IF (IZF-l) 13,14,13 

13 NN=NN-1 
K = K+1 

11 CONTINUE 

14 NP3=NN+3 
N3=NP3 

DO 10 1=1, N3 
CGNV( I ) =0. 

U( I )=0. 

1C V(I)=0. 

B(2)=0.0 
B( 1)=0.0 
C(2)=C.O 
C( 1)=0.0 
D( 2)=0.0 
E(2)=0.0 
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101 



150 

151 



200 

201 



202 

203 

204 
207 



208 
2 50 



251 

252 

253 



254 

255 

256 

2 57 
258 
300 



351 

352 

353 

354 

400 

401 

450 



H(2)=0.0 
DO 101 J = 3,NP3 
H( J)=A( J+K-2) 

CONTINUE 

T=1.0 

SK=10.C=«'*KF 

IF(H(NP3)) 200,151,200 

U(NP3)=0.0 

V(NP3) =0. 

CONV( NP3)=SK 

NP3=NP3-1 

IF(NP3) 51,51,150 

IF(NP3-3) 54,54,201 

PS=0.0 

QS=0.0 

PT=0.0 

QT=0.0 

S=0.0 

REV=1.0 

SK=10.0**KF 

IF(NP3-4) 51,202,203 

R=-H(4)/H(3) 

GO TO 500 
DO 207 J=3,NP3 
IF(H( J) )204,207,204 
S = S+DLOG(DABS(H( J) ) ) 

CONTINUE 
FPN1=N+1 
S=OEXP( S/FPNl i 
DO 208 J=3,NP3 
H( J)=H( J)/S 
CONTINUE 

IF( DABS(H(4)/H( 3 ) ) -DABS ( H ( NP3-1 ) /H( NP3 ) ) ) 250,252, 252 

CONTINUE 

T=-T 

R-(NP3-4)/2 + 3 
DO 251 J=3,M 
S=H( J) 

JJ=NP3-J+3 
H( J)=H( JJ) 

H( JJ)=S 

IF(QS) 253,254,253 
CONTINUE 
P = PS 
Q=QS 

GO TO 300 
HH2=H( NP3-2 ) 

IF(HH2) 256,255,256 
0 = 1.0 
P=-2.0 
GO TO 257 
Q=H(NP3)/HH2 

P=(H(NP3-1 (NP3-3) ) /HH2 

IF(NP3-5)258,550,258 

R=0.0 

CONTINUE 

DO 490 I=1,L 

DO 351 J=3,NP3 

B( J)=H( J)-P’‘B(J-1 )-Q=i'B( J-2) 

C( J) = B( Jl-P’i'CI J-1)-Q*C( J-2) 

CONTINUE 

IF(H(NP3-1) ) 352,400,352 
CONTINUE 

IF (3(NP3-1>) 353,400,353 
AVHB1=DABS(H(NP3-1 ) /B(NP3-1 ) ) 

IF(AVHBl-SK) 450, 354,354 
B(NP3)=H(NP3)-Q*B(NP3-2 ) 

IF(B(NP3) )401, 550,401 
AVHB2=DABS(H(ND3)/5(NP3) ) 

IF(SK-AVHB2) 550,450,450 
DO 451 J=3,NP3 
D( J)=H( J I + R’i^DI J-1 ) 
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456 
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502 



503 
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555 
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E( J)=D( J)+R*E( J-1) 

CONTINUE 

IF(D( NP3) )452 ,500,452 
CONTINUE 

AVH03=DA3S( H(NP3)/D(NP3) ) 

IF ( SK-AVHD3) 500 ,453,453 
CC2=C( NP3-2 ) 

CC3=C( NP3-3 ) 

C(NP3-1 )=-P*CC2-0*CC3 
CC1=C(NP3-1 ) 

S=CC2=!'CC2-CC1*CC3 
IF( 5)455,454,455 
CONTINUE 
P=P-2.00 
Q=Q*(Q+1.0) 

GO TO 456 

P=P+(B (NP3-1 )«CC2-B(NP3)*CC3)/S 
Q=Q+(-B(NP3-l) -CC1+B(NP3)*CC2)/S 
IF(E(NP3-1 ) >458,457,458 
R=R-1.0 
GO TO 490 

R=R-D(NP3)/E(NP3-1 ) 

CONTINUE 
PS=PT 
QS=QT 
PT = P 
QT = Q 

IF(REV >491,492,492 
SK=SK/10.0 
REV=-REV 
GO TO 250 
IF(T) 5C 1 ,502,502 
R=1.0/R 
NP=NP3-3 
U( NP> =R 
NP>-C .0 
CONV(NP>=SK 
NP3=NP3-1 
DO 503 J = 3,NP3 
H(J>=D( J> 

IF(NP3-3>3C0,51,300 

IF(T > 551,552, 552 

P=P/Q 

0=1. 0/Q 

PP2=P/2.C 

0PPSQ=Q-PP2*PP2 

IF( QMPSO >554, 554,553 

NP=MP3-3 

U(NP>=-PP2 

U(NP-1 >=-PP2 

S=DSQRT( CMPSO) 

V(NP>=S 
V(NP-1 >=-S 
GO TO 561 
S=DSQRT(-CMPSQ> 

NP=NP3-3 
IF(P>555,556,556 
U( MP> =-PP2+S 
GO TO 557 
U{NP)=-PP2-S 
U(NP-1>=Q/U(NP> 

V(NP>=0.0 

V( NP-1 > =0.0 

CONV(NP>=SK 

C0NV(NP-1)=SK 

NP3=NP3-2 

DO 558 J=3,NP3 

H( J>=B( J> 

GO TO 200 
IER=1 
RETURN 
END 
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SUBROUTINE ARRAY 



PURPOSE 

CONVERT DATA ARRAY FROM SINGLE TO DOUBLE DIMEN- 
SION OR VICE VERSA. THIS SUBROUTINE IS USED TO 
LINK THE USER PROGRAM WHICH HAS DOUBLE DIMENSION 
ARRAYS AND THE SSP SUBROUTINES WHICH OPERATE 
ON ARRAYS OF DATA IN A VECTOR FASHION. 

USAGE 

CALL ARRAY ( MOO E , I » J t N , M , S , D ) 



DESCRIPTION 

MODE 



I 

J 

N 



D 



=2 - FROM 
NUMBER OF 
NUMBER OF 
NUMBER OF 



OF PARAMETERS 

CODE INDICATING TYPE OF CONVERSION 
=1 - FROM SINGLE TO DUUBLE PRECISION 
DOUBLE TO SINGLE PRECISION 
ROWS IN ACTUAL DATA MATRIX 
COLUMNS IN ACTUAL DATA MATRIX 
ROWS SPECIFIED FOR THE MATRIX 
D IN DIMENSION STATEMENT 
IF MODE=l, THIS VECTOR IS INPUT WHICH 
CONTAINS THE ELEMENTS OF A DATA MATRIX 
OF SIZE I BY J. COLUMN I+l OF DATA 
MATRIX FOLLOWS COLUMN I, ETC. IF MODE=2 
THIS VECTOR IS OUTPUT REPRESENTING A 
DATA MATRIX OF SIZE I BY J CONTAINING 
ITS COLUMNS CONSECUTIVELY. THE LEGNTH 
OF S IS IJ=I*J. 

IF MODE=l, THIS MATRIX OF SIZE N BY M 
IS OUTPUT, CONTAINING A DATA MATRIX OF 
SIZE I BY J IN THE FIRST I ROWS AND J 
COLUMNS. IF MODE=2, THIS N BY M MATRIX 
IS INPUT CONTAINING A DATA MATRIX OF 
SIZE I BY J IN THE FIRST I 
COLUMNS 



ROWS ANU J 



^ ^VECTOR S CAN BE IN THE SAME LOCATION AS MATRIX 
D. VECTOR S IS REFERANCED AS A MATRIX IN OTHER 
SSP ROUTINES, SINCE IT CONTAINS A DATA MATRIX. 
THIS routine converts ONLY GENERAL DATA MATRICES 
( STORAGE MODE 0 ) 



SUBROUTINES 

NONE 



AND FUNCTION SUBPROGRAMS REQUIRED 



SUBROUTINE ARRAYIMODE, I,J,N,M,S,D) 

DIMENSION S(l) ,0(1 ) 

REALMS S,D 

NI=N-I 

TEST TYPE OF CONVERSION 
IF(MODE-l) 100,100,120 

CONVERT FROM SINGLE TO DOUBLE DIMENSION 

100 IJ=I*J+1 
NM=N*J+1 
DO 110 K=1,J 
NM=NM-NI 
DO 110 L = l, I 
IJ=IJ-1 
NM=NM-1 

110 D(NM) = S( IJ ) 
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GO TO 140 

CONVERT FROM DOUBLE TO SINGLE DIMENSION 

120 IJ=0 
NM=0 

DO 130 K=1,J 
DO 125 L=1,I 
IJ=IJ+1 
NM=NM+1 

125 S(IJ)=D(NM) 

130 NM=NM+NI 
140 RETURN 

END 



SUBROUTINE GAUSS3 



PURPOSE 

INVERT A DOUBLE PRECISION MATRIX BY THE GAUSS- 
JORDAN METHOD. THIS ROUTINE IS A DOUBLE PRECIS- 
ION VERSION OF SSP ROUTINE MINV USING Fl-NPGS- 
GAUSS3 (F-63) CALLING SEQUENCE 



USAGE 

DESCRIPTION 

N 

EPS 

A 



KER 



K 



OF PARAMETERS 
ORDER OF MATRIX 
DUMMY PARAMETER 
TV;0 DIMENSIONAL 
TO BE INVERTED 
TWO DIMENSIONAL 
MATRIX 
ERROR FLAG 

=1 INDICATES NO ERRORS 
=2 INDICATES MATRIX IS SINGULAR 
ROW AND COLUMN DIMENSION OF A AND 
USERS PROGRAM 



NOT USED BY GAUSS3 
ARRAY CONTAINING MATRIX 

ARRAY CONTAINING INVERTED 



X IN 



SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

(1) MINV 

(2) ARRAY 

REMARKS 

ALL FLOATING POINT VARIABLES ARE DOUBLE PRECISION 
(REAL^‘8). IF N IS GREATER THAN 50, THE DIMENSION 
OF ARRAYS L,M, AND Y MUST BE CHANGED TO BE 
GREATER THAN OR EQUAL TO N. 



1 



SUBROUT I NE GAUS S3 ( N , EPS , A, X , KER , K ) 

DIMENS lON^M 1) ,X( 1 ) , L ( 50 ) , M ( 50 ) , Y ( 50 , 50 

DO 1 I=1,N 

DO 1 J=1,N 

IND=( I-l )*K+J 

Y(I, J»=A(IND) 

KER=1 



vl2=2*N 

:aLL ARRAY(2,N,N,50,50, Y,Y 

:all minv(y,n,d,l,m) 

:ALL ARRAY! 1 ,N , N, 50 ,50, Y, Y 
IFIO.EQ.O.) KER=2 
DO 2 1=1, N 
DO 2 J=1,N 
IND=( I-l )=«=K+J 



X( IND)=Y( ! ,J) 

RETURN 

END 



) 

) 



) 
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SUBROUTINE MINV 



PURPOSE 

INVERT 

JORDAN 



A DOUBLE 
METHOD. 



PRECISION MATRIX BY THE GAUSS- 



USAGE 

CALL MINV(A,N,D, L,M) 



DESCRIPTION 

A 

N 

D 

L 

M 

subroutines 

NONE 



OF PARAMETERS 

INPUT MATRIX, DESTROYED IN COMPUTATION 

AND REPLACED BY RESULTANT INVERSE. 

ORDER OF MATRIX A 

RESULTANT DETERMINANT 

WORK VECTOR OF LEGNTH N 

WORK VECTOR OF LEGNTH N 

AND FUNCTION SUBPROGRAMS REQUIRED 



REMARKS 

NONE 



C 

C 

c 



c 

c 

c 



SUBROUTINE MI NV ( A, N , D, L , M ) 

DIMENSION A(l) ,L(1) ,M(1) 

DOUBLE PRECISION A , 0 , B I GA , HOLD 

SEARCH FOR LARGEST ELEMENT 

D=1.0 

DO 80 K=1,N 
NK=NK+N 
L( K)=K 
M(K)=K 
KK=NK+K 
B1GA=A( KK) 

DO 20 J=K,N 
IZ=N=^'(J-1 ) 

D0_20 I=K,N 

10 IFIOABS ( BIGA)-DABS( A{ IJ) ) ) 15,20,20 

15 BIGA=A{ IJ ) 

L(K)=I 
M(K) =J 
20 CONTINUE 

INTERCHANGE ROWS 

J=L(K) 

IF(J-K) 35,35,25 
25 KI=K-N 

DO 30 1=1, N 
KI=KI+N 
HOLD=-A(KI ) 

JI=KI^K+J 
A(KI )=A(JI ) 

30 A(JI) =HOLD 

INTERCHANGE COLUMNS 



35 

38 



I=M(K) ^ 

IF(I-K) 45,A5,38 

jP=N*{ i-l ) 

DO 40 J=1,N 

JK=NK+J 

JI=JP+J 
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HOLD=-A( JK) 

A( JK)=A( JI ) 

AO A(JI) =HOLD 

DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT 
ELEMENT IS CONTAINED IN BIGA) 

45 IF(BIGA) 48,46,48 

46 D=0.0 
RETURN 

48 DO 55 1=1, N 

IF(I-K) 50,55, 50 
50 IK=NK+I 

A(IK)=A( IK)/(-BIGA) 

55 CONTINUE 

REDUCE MATRIX 

DO 65 1=1, N 
IK=NK+I 
H0LD=A( IK) 

IJ=I-N 
DO 65 J=1,N 
I J=IJ + N 

IF(I-K) 60,65,60 
60 IF(J-K) 62,65,62 
62 KJ=IJ-I+K 

A( IJ)=HOLD*A(KJ )+A( IJ ) 

65 CONTINUE 

DIVIDE ROW BY PIVOT 

KJ=K-N 

DO 75 J=1,N 

KJ=KJ+N 

IF(J-K) 70-75,70 
70 A(KJ)=A(KJ )/BIGA 
75 CONTINUE 
D=D*BIGA 

REPLACE PIVOT BY RECIPROCAL 

A(KK)=1 .0/8IGA 
80 CONTINUE 

FINAL ROW AND COLUMN INTERCHANGE 



K=N 

100 K=(K-1) 

IF(K) 150,150,105 
105 I=L(K) 

IF(I-K) 120,120,108 
108 JC=N*(K-1) 

JR=N*( I-l ) 

DO 110 J=1,N 
JK=JQ+J 
H0LD=A( JK) 

JI=JR+J 
A( JK)=-A( JI ) 

110 A(JI) =H0LD 
120 J=M(K) 

IF(J-K) 100,100,125 
125 KI=K-N 

DO 130 1=1, N 
K1=KI+N 
HOLD=A( KI ) 

JI=KI-K+J 
A(KI)=-A(JI ) 

130 A(JI) =HCLD 
GO TO ICC 
150 RETURN 
END 
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SUBROUTINE PROD 
PURPOSE 

TO COMPUTE THE PRODUCT OF TWO MATRICES 
USAGE 

CALL PROD(A,B,C,N,M,L) 

DESCRIPTION OF PARAMETERS 

A - FIRST INPUT MATRIX 

B - SECOND INPUT MATRIX 

C - PRODUCT OF A AND B 

N - NUMBER OF ROWS IN A AND C 

M - NUMBER OF COLUMNS IN A AND ROWS IN B 

L - NUMBER OF COLUMNS IN B AND C 

REMARKS 

NONE 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 



SUBROUTINE PROD ( A, B , C , N t M, L ) 
REALMS A(9,9) ,B(9,9) ,C(9,9) 
DO 1 1=1, N 
DO 1 J=1,L 
C( I, J)=0.00 
DO 1 K=1,M 

1 C(I,J)=C(I,J)+A(I,K)>!^B(K,J) 
RETURN 
END 



SUBROUTINE SUMM 
PURPOSE 

ADD TWO MATRICES 
USAGE 

CALL SUMM(A, B,C,M,N) 

DESCRIPTION OF PARAMETERS 

A - NAME OF FIRST INPUT MATRIX 

B - NAME OF SECOND INPUT MATRIX 

C - NAME OF OUTPUT MATRIX 

M -* NUMBER OF ROWS IN A,B, AND C 

N - NUMBER OF COLUMNS IN A,B, AND C 

REMARKS 

NONE 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 



SUBROUTINE SUMM ( A , B , C , M , N ) 
REALMS A(9,9) ,B(9,9),C(9,9) 
DO 1 1=1, M 
DO 1 J=1,N 

1 C(I,J) = A(I, J)+B( I, J) 

RETURN 

END 
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SUBROUTINE DIFF 
PURPOSE 

SUBTRACT ONE MATRIX FROM ANOTHER 
USAGE 

CALL DIFF(A, B,C,M,N) 

DESCRIPTION OF PARAMETERS 

A - FIRST INPUT MATRIX 

B - SECOND INPUT MATRIX 

C - OUTPUT MATRIX EQUALS A - B 

M - NUMBER OF ROWS IN A,B, AND C 

N - NUMBER OF COLUMNS IN A,B, AND C 

REMARKS 

NONE 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 



SUBROUTINE D IF F ( A, B t C, M , N ) 
REALMS A(9,9) ,B(9,9 ) ,C(9,9) 
DO 1 I=1,M 
DO 1 J=1,N 

1 C(I,J)=A(I, J)-B(I, J) 

RETURN 

END 



SUBROUTINE EXPAND 
PURPOSE 

TO COMPUTE THE REAL COEFFICIENTS OF AN N-TH 
DEGREE POLYNOMIAL GIVEN N COMPLEX ROOTS 

USAGE 

CALL EXPAND(N,R, A) 

DESCRIPTION OF PARAMETERS 

N - DEGREE OF POLYNOMIAL 

R - VECTOR CF COMPLEX ROOTS 

A - COEFFICIENT VECTOR 

REMARKS 

NONE 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 



1 

2 

3 



SUBROUTINE EXPAND! N, R, A) 
REALMS A(10) 

CCMPLEX=i=16 R(9 ) ,Q(9) 

IF (N-1) 1,2,3 

A(l)=1.00 

RETURN 

A( 1 )=-REAL(R(l) ) 
A(2)=1.00 
RETURN 
Q(l)=-R(l) 

Q(2)=1.00 
DO 5 J=2,N 
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Q( J+l)=1.00 
JJ=J-1 
DO 4 I=1»JJ 
K=J-I 

4 Q(K+1)=Q(K)-Q(K+1 )«R( J) 

5 Q(1)=-Q(1)=!'R(J) 

6 N1=N+1 

DO 7 L=1,N1 

7 A(L)=REAL(Q(L) ) 

RETURN 

END 
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